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ABSTRACT 

The  processing  of  thermoset  energetic  polymers  involves  processes  in  which  a 
polymerization  of  an  initial  liquid  prepolymer  takes  place.  This  chemical  reaction  results 
in  a  significant  increase  of  the  molecular  weight  of  the  polymer.  Consequently,  the 
rheological  properties  of  the  material  are  also  affected  by  this  curing  phenomenon.  The 
extent  of  this  reaction,  as  well  as  the  type  of  polymer  used,  often  determine  the 
magnitude  of  these  changes  and  their  effect  on  the  quality  of  the  end-product  being 
processed.  In  order  to  optimize  the  processing  window  of  polyurethane-based 
formulations,  the  polymerization  reaction  has  been  modeled  by  Monte-Carlo  simulations 
For  this  purpose,  a  numerical  code  has  been  developed  in  ANSI  Fortran  77  that  allows 
the  simulation  of  A2+B2  polymerization  with  provision  for  unequal  reactivity  of  the 
reaction  sites.  Simulations  have  been  carried  out  on  HTPB-TDI,  PPG-HDI  and  GAP- 
IPDI  systems.  A  limited  experimental  validation  has  confirmed  the  validity  of  the 
molecular  weight  distributions  calculated  by  the  numerical  code.  The  information  obtained 
from  these  simulations  should  provide  a  means  to  estimate  the  material  functions  of  the 
polymer  during  the  mixing  and  casting  of  a  formulation  by  transforming  the  molecular 
weight  distribution  into  a  relaxation  time  distribution.  Minor  modifications  to  the 
algorithm  will  enable  the  simulation  of  more  complex  An+Bn  systems. 


RESUM& 


La  mise  en  oeuvre  de  polymdres  dnergdtiques  thermodurcissables  implique  des 
operations  comprenant  la  polymerisation  d’un  prepolymere  liquide.  Cette  reaction 
chimique  se  traduit  par  une  augmentation  importante  de  la  masse  molaire  du  polymdre.  En 
consequence,  les  propridtds  rhdologiques  du  materiel  sont  elles  aussi  trds  affectdes  par  le 
phdnomdne  de  cuisson.  L'avancement  de  la  reaction,  ainsi  que  le  type  de  polymdre  utilise, 
determinent  souvent  l'importance  de  ces  changements  et  leur  effet  sur  la  qualitd  des 
produits  finis.  Afin  d’optimiser  les  conditions  d’operations  pour  la  mise  en  oeuvre  de 
compositions  k  base  de  polyurdthanne,  la  reaction  de  polymerisation  a  did  moddlisde  par 
des  simulations  numdriques  stochastiques.  A  cet  effet,  un  code  a  dtd  ddveloppd  en  Fortran 
77  selon  la  norme  ANSI.  Ce  programme  permet  la  simulation  de  systdmes  de  type 
A2+B2  avec  la  possibilitd  de  prendre  en  consideration  la  difference  de  rdactivitd  entre  les 
sites.  Des  simulations  ont  dtd  effectudes  pour  les  systdmes  rdactifs  HTPB-TDI,  PPG- 
HDI  et  GAP-IPDI.  Des  travaux  expdrimentaux  ont  permis  de  valider  les  rdsultats  de 
distribution  des  masses  molaires  obtenus  par  le  logiciel  de  simulation.  II  appert  que 
l’information  obtenue  permettrait  d’estimer  les  fonctions  matdrielles  du  polymdre  durant  le 
malaxage  et  la  coulde  d’une  composition  en  transformant  la  distribution  des  masses 
molaires  en  distribution  des  temps  de  relaxation.  Des  modifications  mineures  devront  etre 
apportdes  au  code  afin  de  permettre  la  simulation  de  systdmes  plus  complexes  de  type 
An+Bn- 


P154742.PDF  [Page:  6  of  78] 


P154742.PDF  [Page:  7  of  78] 


UNCLASSIFIED 

iii 

aRECEDiHG  PAGE  seank 

TABLE  OF  CONTENTS 


abstract/resum£ .  i 

EXECUTIVE  SUMMARY .  v 

LIST  OF  SYMBOLS .  vii 

1.0  INTRODUCTION .  1 

2.0  THEORY 

2. 1  Description  of  the  Time  Behavior  of  a  Chemical  Reaction .  2 

2.2  Effect  of  Molecular  Weight  on  Rheological  Properties .  3 

2.3  Statistical  Estimation  of  Lower  Moments  of  the  Molecular 

Weight  Distribution  .  6 

2.4  Stochastic  Algorithm  Describing  Polymerization  Reaction  .  9 

3.0  CODEDEVELOPEMENT 

3.1  Molecular  Digitization .  11 

3.2  Performing  a  Reaction .  14 

3.3  Algorithm  and  Software .  15 

4.0  SIMULATIONS 

4.1  Comparative  Polyurethane  Systems .  17 

4.2  Effect  of  Polymerization  Reaction  Variables .  19 

5.0  EXPERIMENTAL  .  22 

6.0  DISCUSSION .  23 

7.0  CONCLUSIONS .  27 

8.0  ACKNOWLEDGEMENTS .  28 

9.0  REFERENCES . . . 29 

TABLES  I  to  IV 
FIGURES  1  to  12 

APPENDIX  A-  Listing  of  the  Computer  Code  for  MWD  Estimations 
APPENDIX  B-  Template  for  Input  Data  Files 
APPENDIX  C-  Example  of  Screen  Output 


P154742.PDF  [Page:  8  of  78] 


P154742.PDF  [Page:  9  of  78] 


UNCLASSIFIED 

fiRECEOJNG  PAGE  SiANK 

EXECUTIVE  SUMMARY 

In  times  of  constant  reduction  of  spending  for  defence  purposes,  military 
industries  are  challenged  to  produce  state-of-the-art  weapons  and  ammunition  at  the  lowest 
cost  possible.  Hence,  any  technology  transferred  from  CRAD  to  industry  must  be 
developed  in  keeping  in  mind  these  economic  considerations.  Accordingly,  research  and 
development  in  explosives  and  propellants  must  be  carried  out  in  such  a  way  that  new 
formulations  can  be  easily  produced  by  Canadian  industry.  For  thermoset  composite 
materials,  the  knowledge  of  the  rheology  of  a  formulation  is  critical  information  in  the 
establishment  of  a  processing  window  for  mixing  and  casting  operations. 

Recent  advances  in  polymer  rheology  have  demonstrated  the  possibility  of 
estimating  most  of  the  material  functions  from  the  molecular  weight  distribution  (MWD)  of 
the  polymer.  Unfortunately,  the  molecular  weight  distribution  of  thermoset  polymers 
changes  continuously  over  the  processing  cycle,  starting  from  the  moment  that  the  curing 
agent  is  added  to  the  formulation.  In  a  first  step  toward  a  predictive  model  for  explosive 
and  propellant  production  processes,  this  memorandum  reports  the  use  of  a  stochastic 
algorithm  to  describe  quantitatively  the  evolution  of  molecular  species  during  polyurethane 
formation.  A  computer  code  was  developed  to  carry  out  this  algorithm.  Simulations  on 
simplified  formulations  based  on  hydroxy-terminated  polybutadiene  (HTPB)  and  glycidyl 
azide  polymer  (GAP)  have  been  performed.  The  molecular  weight  distributions  obtained 
have  been  verified  by  gel  permeation  chromatography  experiments.  From  these  results,  it 
appears  that  the  software  would  be  able  to  predict  the  molecular  weight  distribution  of  most 
polyurethane  formulations  before  gel  point. 

Since  a  significant  linkage  exists  between  the  rheology  of  the  polymer  and  its 
predicted  molecular  structure,  the  Monte-Carlo  simulations  would  be  a  very  useful  tool  for 
the  development  of  new  propellant  formulations.  By  simulating  the  processing  of  new 
thermoset  polymers,  the  cost  of  their  development  would  be  reduced  and  more  ingredient 
combinations  could  be  tested  in  the  same  amount  of  time.  Also,  formulations  that  could  not 
be  processed  into  a  defect-free  end-product  would  be  discarded  from  the  start.  Hence,  more 
complex  energetic  materials  will  be  developped  and  this  will  benefit  to  the  Canadian  Armed 
Forces  by  making  available  to  them  more  powerful  weapon  systems,  at  a  better  cost 
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1.0  INTRODUCTION 

The  increasing  use  of  composite  materials  in  industries  today  continuously  raises 
new  challenges  for  the  processing  of  thermoset  polymers.  In  the  field  of  military 
technology,  composite  products  are  found,  besides  the  well-known  aerospace 
applications,  in  many  energetic  materials.  These  include  polyurethane-bonded  composite 
explosives  and  solid  rocket  propellants.  The  development  of  new  propellant  formulations 
requires  a  thorough  evaluation  of  all  aspects  of  the  material  characteristics.  One  of  them 
concerns  the  ability  of  the  propellant  to  be  casted  into  a  defect-free  propellant  grain.  In 
order  to  achieve  this  requirement,  special  care  must  be  taken  to  optimize  the  processing 
cycle.  When  casting,  the  paste  must  retain  enough  fluidity  to  prevent  void  formation, 
while  showing  a  sufficient  consistency  to  avoid  sedimentation. 

Among  the  material  functions  of  interest  in  the  evaluation  of  solid  rocket 
propellants,  the  viscosity-shear  rate  spectrum  is  of  utmost  importance  (Ref.  1). 
Numerous  techniques  have  been  studied  to  correlate  the  observed  rheology  of  the  curing 
polymers  with  the  progress  of  the  polymerization  reaction.  Often,  this  chemorheological 
characterization  is  performed  by  adiabatic  viscometry  (Ref.  2).  In  such  a  case,  the 
change  in  rheology  is  monitored  and  correlated  with  time.  Unfortunately,  these  relations 
cannot  be  easily  adaptated  for  changes  of  the  polymerization  process  parameters,  such  as 
temperature  or  catalyst  concentration.  In  some  instances,  the  viscosity  has  been  reported 
as  a  function  of  the  extent  of  reaction  (Refs.3,4).  This  approach  makes  possible,  with 
appropriate  kinetics  expressions,  the  prediction  of  how  changes  in  the  reaction  conditions 
will  effect  the  conversion  and,  therefore,  the  rheology  of  the  polyurethane  system. 

However,  the  kinetic  (or  the  deterministic)  approach  alone  is  unable  to  self¬ 
adjust  to  changes  in  the  prepolymer  to  be  cured.  So  different  polyols  will,  for  the  same 
level  of  conversion,  obviously  have  different  rheologies,  simply  because  their  molecular 
weights  will  be  different  This  leads  to  the  uneasy  task  of  estimating  the  molecular  weight 
of  the  polymer  as  a  function  of  the  extent  of  reaction.  In  fact,  the  whole  molecular  weight 
distribution  must  be  found  if  one  wants  to  obtain  truly  complete  rheological  information. 
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In  this  memorandum,  an  application  of  stochastic  simulation  algorithms  to  the 
reaction  of  polyurethane  systems  found  in  rocket  propellants  is  presented.  The  objective 
of  the  work  described  here  is  to  demonstrate  the  usefulness  of  this  methodology  as  a 
means  of  modeling  the  processing  of  reactive  polymeric  material.  In  this  present  study, 
the  scope  of  the  modeling  is  restricted  to  conversion  up  to  the  gel  point.  A  numenc  code 
in  Fortran77  has  been  developed  that  allows  calculations  on  multi-component  polyol- 
polyisocyanate  mixtures.  The  core  of  the  software  will  run  on  any  platform  being  able  to 
compile  ANSI  FORTRAN77  code,  while  the  graphics  interface  modules  included  are 
specific  to  Apple  Macintosh®  operating  systems.  Experiments  have  been  carried  out  to 
validate  the  numeric  simulations.  Three  different  were  systems  considered:  HTPB-TDI, 
PPG-HDI  and  GAP-IPDI. 

This  technical  memorandum  reports  work  done  at  DREV  between  January  1994 
and  December  1994  under  project  PSC  32C,  Missile  Propulsion. 


2.0  THEORY 

2.1  Description  of  the  time  behavior  of  a  chemical  rgact19.il 


The  problem  of  describing  the  evolution  of  a  chemically-reacting  system  is 
mathematically  addressed  by  establishing  appropriate  time  functions  for  concentrations  of 
each  species  found  in  the  reaction  space.  As  Gillespie  (Ref.  5)  pointed  out,  this  is 
generally  accomplished  through  either  a  deterministic  or  a  stochastic  approach. 

The  deterministic  point  of  view  is  the  more  usual  way  of  describing  the  kinetics  of 
a  chemical  reaction.  Generally,  the  rate  of  change  for  the  population  of  each  of  the  n 
components  of  a  system  is  expressed  by  an  equation  like: 

. x.)  111 

at 

Hence,  the  overall  behavior  of  the  reaction  is  given  by  a  set  of  n  ordinary  differential 
equations.  Often,  an  analytical  solution  of  the  problem  is  impossible  to  achieve  and 
extensive  numerical  methods  have  to  be  used  (Ref.  6)  to  solve  equation  [1]. 
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The  effectiveness  of  the  reaction  rate  equations  to  evaluate  the  kinetics  has  been 
demonstrated  for  an  impressive  number  of  reactive  systems.  However,  it  must  be  kept  in 
mind  that,  mathematically,  this  method  interpolates  a  continuous  function  over  a 
physically-discrete  phenomenon.  Indeed,  the  molecular  populations  of  a  reaction  space 
cannot  evolve  by  a  quantity  other  than  the  number  (integer)  of  molecules  taking  part  in  a 
reaction.  Fortunately  though,  for  most  problems,  this  quantum  limitation  can  be 
overlooked. 

The  stochastic  approach  considers  the  evolution  of  a  chemical  reaction  as  the 
macroscopic  result  of  a  purely  random  series  of  events  taking  place  at  the  molecular  scale. 
Fundamentally,  the  description  of  the  time-dependent  process  is  given  by  a  distribution  of 
probability  functions  for  each  elementary  reaction  occurring  within  the  reaction  space. 
These  functions  will  have  to  be  defined  by  simple  formulations  based  on  the  levels  of 
species  forming  the  reaction  space  population.  Whenever  possible,  a  prior  knowledge  of 
the  reaction  pathway  will  be  helpful.  The  probability  distributions  are  therafter  used  either 
in  a  master  probability  equation  or  in  Monte-Carlo  simulations.  This  shall  be  discussed  in' 
section  2.4.  Examples  of  the  use  of  the  stochastic  determination  of  reaction  kinetics  have 
been  reported  to  solve  problems  related  to  inorganic,  pharmaceutical  and  polymer 
materials  chemistry  (Refs.  5-10). 

2.2  Effect  of  molecular  weight  on  rheological  properties 

There  is  an  intimate  relation  between  the  size  of  polymer  chains  in  the  liquid  state 
(or  concentrated  solution)  undergoing  shear  flow  and  their  material  functions.  The  most- 
often  studied  of  these  rheological  properties,  viscosity,  is  affected  by  the  rate  of  shear 
applied  on  the  material.  Typically,  for  the  range  of  shear  rate  found  in  polymer  mixing 
processes,  the  viscosity  -  shear  rate  spectrum  plotted  on  a  logarithmic  scale  will  have  a 
plateau  portion  at  lower  shear  rates  followed  by  a  linear  negative  slope.  These  two 
sections  of  the  spectrum  are  often  referred  to  as  the  zero  shear  rate  viscosity  (Newtonian) 
and  the  power-law  index  (shear  thinning).  Both  are  affected  by  the  size  of  the  polymer 
molecule.  The  zero  shear  rate  viscosity  (r\0)  changes  as  a  power  of  the  molecular  weight 
(M).  For  monodisperse  polymers,  the  following  relation  holds  (Ref.  1 1): 
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rio~Ma 


[2] 


where  1  <  a  <  2.5  for  M  <  Me 
a  =  3.4  for  M  >  Me 


It  is  seen  that  beyond  the  molecular  weight  of  entanglement  (Me)  the  relation  between  rio 
and  M  will  always  be  the  same,  regardless  of  the  considered  species.  Accordingly,  a  fully 
logarithmic  plot  of  the  zero  shear  rate  viscosity  against  the  molecular  weight  will  exhibit 
two  straight  lines  with  different  slopes  crossing  at  Me.  For  polydisperse  polymers,  the 
weight-average  molecular  weight  (Mw)  of  the  species  shall  be  used  in  equation  [2]  (Ref. 
12).  These  relations  have  proved  to  be  remarkably  suitable  for  a  very  large  number  of 
linear  or  slightly  branched  polymers. 


Even  though  the  zero  shear  rate  viscosity  is  the  primary  material  function  used  for 
characterization  and  simulation  purposes,  the  dependence  of  t\q  on  shear  rate  must  also  be 
known  if  one  wants  to  describe  a  chemical  process  where  a  wide  spectrum  of  shear  rates 
are  applied  to  the  polymer,  as  with  rocket  propellant  mixing.  In  such  a  case,  it  may  be 
useful  to  describe  the  transition  between  the  Newtonian  regime  to  the  shear  thinning 
behavior  by  a  characteristic  time  (Kq)  based  on  molecular  weight  ( Refs.  12-14): 


A  =TJ°M 
e  pRT 


[3] 


Again,  unfortunately,  some  limitations  arise  from  the  fact  that  the  slope  of  the  power-law 
region  may  change  with  Ae,  in  such  a  way  that  a  theoretical-based  correlation  is  difficult 
to  establish. 


Recently,  Tsenouglou  (Ref.  15)  and  Mead  (Refs.  16,  17)  have  studied  the 
relationship  between  molecular  weight  and  relaxation  time  spectrum  of  polydisperse 
homopolymers.  Mead  has  shown  that  equivalent  information  was  stored  in  the  molecular 
weight  distribution  (MWD)  of  a  polymer  and  in  its  time  relaxation  spectrum,  given  that 
numerical  interconversion  was  possible.  This  can  be  accomplished  by  the  following 
procedure: 
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1-  simulate  (or  assume,  or  measure)  a  MWD 

2-  evaluate  the  relaxation  time  spectrum  (  GbUnd  ( t ) )  by: 

=  J~  W(M)4¥(Mj)dM  [4] 

where  *jF(M,t)  is  the  monodisperse  species  relaxation  time  function.  Often,  when  no 
experimental  data  are  available,  this  is  expressed  by  an  exponential  of  the  form: 


s]F(M,t)  [5] 

3-  from  the  relaxation  spectrum,  the  dynamic  material  functions  are  easily  found  by: 

G  (co)  =  mj;  Gblend  (t)sm(cot)dt  [6] 

G"  (co)  =  co  j~  Gblend  (i t )  cos(cot)dt  [7] 

G(co)  =  i}"co  [8] 

G"(co)  =  r[co  [9] 


4-  finally  the  steady  shear  viscosity  is  estimated  by  the  Cox-Mayer  rule  (Ref.  12)  and  the 
definition  of  the  zero  shear  rate  viscosity: 


77=77*  =  77'  (co) 


■ 

(  ntf  ^ 

2" 

1+ 

U’J 

-10.5 


N 


77  o=  Jo°  G(t)  dt  =51  Gfii 


[10] 

[11] 


Of  course,  the  completion  of  this  procedure  requires  numerical  methods  for  the  solving  of 
Fredholm  integral  equations  of  the  first  kind  that  will  not  be  presented  into  this 
memorandum.  However,  it  shall  be  seen  that  the  primary  motivation  for  evaluating  the 
MWD  through  Monte-Carlo  simulation  lays  in  the  potential  rheological  information 
generated  by  this  process. 
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It  has  been  established  that  a  function  like  a  molecular  weight  distribution  can  be 
rigorously  described  by  all  of  its  integration  moments  (Ref.  17).  However,  for  polymer 
science  and  polymer  processing,  only  the  fourth  first  moments  are  usually  retained  for 
characterization  purposes.  These  are: 

i)  zero  moment: 

Mo=in;M;  [12] 

i 

it  simply  states  that  the  molecular  weight  distribution  must  be  closed  for  the  species' 
mass  fractions. 


ii)  first  moment: 

N 

2>«- Mi 

Mn^-hr— 

Ini 


often  referred  to  as  the  number-average  molecular  weight 


iii)  second  moment: 

N 

ln,M? 

tniMi 

i 


usually  known  as  the  weight-average  molecular  weight 


[13] 


[14] 


iv)  third  moment: 
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N 

ZriiM? 

-  [15] 

i 

simply  called  the  Z-  average  molecular  weight  (Mz) 

The  determination  of  these  moments  of  the  MWD  as  a  function  of  the  advancement  of  a 
polymerization  reaction  has  been  studied  via  different  methods  over  the  last  40  years. 
Some  of  them,  found  in  open  literature,  will  be  presented  below. 

The  zero  moment  is  clearly  trivial  to  verify.  The  first  moment  can  always  be 
calculated  by  simple  stoechiometry.  For  example,  one  may  consider  a  polymerization 
reaction  where: 


A2  +B2  -->  high  molecular  weight  polymer 


where  A2  and  B2  are  a  diol  prepolymer  and  a  difunctional  isocyanate  respectively.  This 
results  in  the  formation  of  a  polyurethane  .  The  number-average  molecular  weight  at  a 
particular  instant  of  the  reaction  is  the  total  mass  over  the  total  number  of  molecules.  For 
the  previous  reaction  this  translates  by: 


where: 


and 


Ma  +  Mn  r 

M  = _ -l _  2 

"  l+r-2X 


nB2 


X  = 


„0 


[16] 


[17] 


[18] 


The  second  moment  of  MWD  requires  a  more  challenging  analysis  for  its  determination. 
Early  researchers  such  as  Stockmayer  (Ref.  18)  used  extensive  algebraic  manipulations 
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based  on  probability  distributions  to  deduce  a  typical  expression  for  polycondensation. 
Keeping  the  same  example  for  simplification,  the  appropriate  expression  is: 


(1 Jr) 


Mw=- 


1  + 


X 


2\ 


*1  + 


1  +  - 


.  4 XMAMB2 


m2  + 


rM%  0 


V 


1- 


X 


2\ 


[19] 


A 


In  the  case  of  a  polydisperse  prepolymer,  Ziegel  et  al  (Ref.  19)  found  that  this  relation 
can  be  expressed  by: 


m 


Mw=- 


i- 


X 


2  'N 


Mn  M  + 
nA2  wa2 


(  y2  ^ 

1-  — 
r 


Ml  +— (xM^  +~Ml  +2M„ .  Mr, 


^~nA2 


lnA2lViB2 


m ; 


2-  + M 


Bn 


/l_ X*^ 


J 


V 


J 


[20] 

Later  Macosko  et  al  (Refs.  20,  and  21)  carried  out  a  generalization  of  this  procedure  with 
a  much  simpler  analysis  based  on  the  expectancy  of  probabilities  for  a  defined 
configuration. 

The  moments  predicted  by  these  mathematical  techniques  have  been  confirmed  for 
many  polymeric  systems  by  experimental  analytical  tools  like  GPC,  tonometry, 
viscometry,  etc.  They  are  particularly  useful  for  estimation  of  gel  points  in  polymerization 
reactions.  Unfortunately,  they  are  not  as  reliable  for  high  degree  of  conversion.  In  such 
situations,  the  mechanisms  of  molecular  species  diffusion  must  be  considered.  Gupta  et  al 
(Ref.  10)  have  demonstrated  the  use  of  Monte-Carlo  simulations  in  polymerization 
reactions  where  diffusion  must  be  accounted  for.  In  these  simulations,  the  spatial 
arrangement  of  each  species  was  stored  in  memory  in  order  to  assess  the  probability  of  a 
reaction  between  two  reactive  sites  based  on  the  distance  between  them.  Unfortunately 
this  requires  a  fairly  large  amount  of  computational  effort  The  approach  presented  in  this 
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memorandum  is  similar  with  regards  of  the  Monte-Carlo  algorithm  but  the  computation 
task  is  greatly  reduced  by  the  equireactivity  assumption  since  no  records  of  the  molecule 
positions  are  kept.  However,  time  progress  of  the  reaction  is  monitored. 

2.4  Stochastic  algorithm  describing  polymerization  reaction 

The  simulation  approach  presented  here  makes  use  of  a  Monte-Carlo  technique.  In 
order  to  achieve  a  description  of  a  whole  polymerization  process,  elementary  steps,  or 
events  defining  it  have  to  be  known.  After  that,  the  relative  probability  of  each  event  has 
to  be  evaluated.  Successive  repetition  of  randomly  chosen  events  suffice  to  describe  the 
polymerization  process. 

Consider  the  reaction  described  in  section  2.3,  where  one  of  the  two  functional 
groups  of  the  isocyanate  has  a  different  reactivity.  In  such  a  case,  three  events  can  take 
place: 

-  a  reaction  between  an  OH  group  and  a  NCO  group  (event  1) 

-  a  reaction  between  an  OH  group  and  a  NCO’  group  (event  2) 

-  no  reaction  at  all;  only  time  elapses  (event  3). 

Accordingly,  each  event  will  have  to  be  selected  upon  its  probability,  or  frequency,  of 
occurrence.  If  a  random  number  is  distributed  between  0  and  1,  then  three  intervals  can 
be  defined  as: 

[0,pl]  ,  [pl,pl+p2]  and  [pl+p2,  1] 


where 


Pi  = 


_ P 'event  1 _ 

P event  1  P event  2  P event  3 


[21] 


and  similarly  for  p2  and  p3.  It  results  that  if  R,  a  random  number  describing  the  ith 
event,  is  part  of  the  interval  [0,pl],  then  event  1  will  occur,  and  so  on.  The  difficult  task 
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has  to  go  back  to  traditional  kinetics  information.  Hence,  the  apparent  rate  constants  can 


make  very  acceptable  approximations  of  event  probabilities.  For  the  situation  described 
before,  the  apparent  reaction  rate  equations  are  expressed  by  (Ref.  22) : 

^21l=*1[oh][nco]1 

[22] 

[23] 

Considering  that  a  single  reaction  event  involves  1/F  mole  of  each  species  within  a 

reaction  space  of  volume  V,  then 

A[NCO],  =  A[NCO]2 

[24] 

and  the  incremental  time  for  each  event  is  given  by: 

1 

tl  ife1FVr[OH][NCO]1 

[25] 

-  1 
h  *,FV[OH][NCO]2 

[26] 

At,  -  arbitrary  constant 

[27] 

It  is  obvious  that  the  most  probable  event  should  have  the  smaller  time  increment,  so  the 
normalized  probability  of  events  for  this  reaction  scheme  will  be  expressed  by: 

1/A  t. 

Pi  ~  l/A/j  +  1/At,  +  l/Afj 
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[28] 
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3.0  CODE  DEVELOPMENT 


3.1  Molecular  digitization 

The  simulation  of  a  polymer  chemical  reaction  by  stochastic  modeling  asks  for  a 
translation,  or  an  encoding,  of  the  molecular  structure  information  to  a  numerical  data 
bank  to  be  manipulated  by  the  code.  For  this  work,  an  approach  adapted  from  Pandit  et 
al.  employing  a  well-known  methodology  (Ref.  7  )  has  been  retained.  It  consists  of 
expressing  the  material  domain  by  three  vectors:  the  reaction  space  vector  (RSV),  the 
species  weight  vector  (SWV),  and  the  link  weight  vector  (LWV).  These  quantities  are 
described  below. 

3.1.1  The.jeaelioD.spag,g..ve.g.t.Qr 

The  reaction  space  vector  contains  the  numerical  translation  of  all  the  information 
about  the  molecular  structure  of  every  species  found  in  the  reacting  system.  This  is  done 
by  assigning  a  number  to  each  kind  of  elementary  group  forming  a  molecular  species.  An 
elementary  group  is  considered  to  have  no  internal  links  affected  by  a  chemical  reaction. 
For  example,  an  oligomer  of  polypropylene  glycol  (PPG)  like  : 

H0-CH(CH3)(CH2)0-CH(CH3)(CH2)0-CH(CH3)(CH2)0-CH(CH3)(CH2)-0H 

would  have  a  numeric  equivalent  of 

{0,4  x  4,1}  =  {0,16,1} 

.  The  conversion  from  molecular  to  vectorial  representation  is  performed  by  using  the  key 
described  in  Table  I.  It  is  seen  that  each  hydroxyl  group  is  translated  to  the  numeric  value 
of  0  or  1,  while  successive  monomer  body  groups  expressed  by  4's  are  summed  to 
reduce  the  length  of  the  RSV. 

Inversely,  a  form  like  {0,12,3,5,2}  would  correspond,  for  a  PPG-HDI  system,  to  a 
molecule  like 


H0-CH(CH3)(CH2)0-CH(CH3)(CH2)0-CH(CH3)(CH2)0-NH-C0-0-(CH2)6-NC0 
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Finally,  the  RSV  will  be  constructed  by  sequentially  adding  each  species  vector. 
Henceforth,  for  a  reaction  space  containing  at  some  instant  of  the  polymerization  process 
these  five  chemical  species: 

{0,12,1},  {1,12,3,5,2},  {2,5,2},  {2,5,3,16,3,5,2}  and  {1,8,3,5,3,16,0} 
the  resulting  RSV  is: 

{0,12,1,  1,12,3,5,2,  2,5,2,  2,5,3,16,3,5,2,  1,8,3,5,3,16,0} 

It  shall  be  noted  that  any  particular  species  of  the  reaction  space  vector  can  be  tracked 
back  by  using  the  fact  that  each  of  them  must  end  by  an  OH  or  a  NCO  group. 

TABLE  I 

Numerical  Conversion  Key  lor  jypic&LSy  stems 


Group  Number 

Molecular  Group 
for  HTPB-TDI 
system 

Molecular  Group 
for  PPG-HDI 
system 

Description 

0 

- 

OH 

1 

OH 

OH 

2 

NCO 

3 

NH-CO-O 

NH-CO-O 

urethane  | 

4 

CH2CHCHCH2 

CH(CH3)(CH2)0 

HSBBB 

5 

(CH2)6 

6 

NCO  ortho 

- 

BBS5SS9BI 

3.1.2  The  species  weight  vector 

The  species  weight  vector  contains  the  amount  (in  moles)  of  each  species  found  in 
the  RSV,  keeping  the  same  relative  position  in  the  vectors.  Therefore  the  third  element  of 
the  SWV  is  the  mass  of  the  third  species  found  in  the  RSV. 
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A  typical  SWV  for  the  RSV  expressed  before  would  be: 

{  1.1,  0.09,  1.0,  0.12,  0.3} 

3.1.3  The  link  weight  vector 

The  link  weight  vector  completes  the  set  of  reactor  space  descriptors.  It  consists 
of  the  molar  amount  of  each  particular  chemical  link  between  elementary  groups  used  to 
build  the  RSV.  Accordingly,  if  five  different  kinds  of  links  are  found,  five  elements  will 
constitute  the  LWV.  Table  n  summarizes  chemical  links  to  be  monitored  (numerically) 
for  two-component  polyurethane  systems.  The  LWV  is  calculated  directly  from  the 
information  encoded  into  the  RSV  and  the  SWV.  For  the  previous  example,  and  using  the 
key  of  Table  n,  the  resulting  LWV  is 

{1.4,  1.49,  4.24,  2.3,  0.0,  0.93,  0.93} 


TABLE...P 

Equivalence  Table  for  Chemical  Linkage 


Position 
in  the  LWV 

Groups  Involved 

Link  type 

1 

0-4 

prim,  hydroxyl  -  monomer  body 

2 

1-4 

sec.  hydroxyl  -  monomer  body 

3 

4-4 

monomer  -  monomer  body 

4 

2-5 

isocyanate  -  isocyanate  body 

5 

5-6 

isocyanate  body  -  isocyanate 

6 

3-4 

urethane  -  monomer  body 

7 

3-5 

urethane  -  isocyanate  body 
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3.2  Performing  a  reaction 

Within  the  frame  of  this  particular  study,  only  condensation  reactions  from  an 
hydroxyl  group  with  an  isocyanate  function,  to  form  an  urethane  linkage,  have  been 
considered.  Consequently,  the  total  number  of  reactions  of  different  kinds  susceptible  to 
take  place  into  the  reaction  space  is  directly  related  to  the  number  of  different  (non- 
chemically  equivalent)  OH  and  NCO  functions  among  all  the  species.  For  a  system  with 
Noh  different  hydroxyl  functions  and  Nnco  different  NCO  functions  the  total  number  of 
chemical  reactions  is  simply: 

Ncr  =  Noh  Nnco  [29] 

Often,  hydroxyl  functions  of  a  polydisperse  polymer  are  considered  to  be  equivalent, 
with  no  regards  to  the  molecular  weight  of  a  particular  species  of  the  molecular  weight 
distribution  (MWD).  In  such  a  case,  the  number  of  possible  chemical  reactions  will  be 
imposed  by  the  isocyanate  population.  For  example,  with  the  HTPB-TDI  system,  the  TDI 
species  bears  two  different  NCO  groups,  since  the  reactivity  of  the  ortho  group  is  much 
less  than  that  of  the  para-  position  NCO  function.  Hence,  there  are  two  possible  reactions 
to  consider.  It  has  been  discussed  in  section  2.4  how  the  choice  of  the  reaction  occurring 
at  a  step  of  the  simulation  path  is  performed. 

Once  the  reaction  type  has  been  selected,  the  reaction  sites  have  to  be  chosen.  This 
is  done  on  the  basis  that  the  probability  of  having  a  species  reacting  with  another  is 
directly  proportional  to  the  amount  (mass)  of  this  species  found  into  the  reaction  space. 
So,  if  a  OH  -  NCO  reaction  is  to  be  performed,  it  is  likely  that  the  molecular  species  of 
the  polydisperse  polymer  having  the  highest  mass  fraction  has  more  chance  to  react.  But 
this  does  not  mean  that  it  will  actually  be  the  selected  site  for  the  reaction  (since  the 
probability  of  having  a  particular  species  involved  in  a  specific  reaction  is  always  less 
than  one). 

After  the  selection  of  the  two  reacting  species,  a  copy  of  their  molecular  vector 
(MV)  is  extracted  from  the  RSV: 


Molecule  bearing  the  hydroxyl  reaction  site: 


{1,12,1} 
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Molecule  bearing  the  isocyanate  reaction  site:  {2,5,2.} 

The  chemical  reaction  between  these  sites  is  numerically  simulated  by  forming  a  new 
molecular  vector  incorporating  the  new  link  formed: 

New  molecule:  {2,  5,  2,  12,  1} 

If  the  molecule  formed  is  new  to  the  reaction  space,  it  will  be  added  in  the  RSV  at  an 
index  position  based  on  its  molecular  weight.  After  that,  the  elements  of  the  SWV  and 
LWV  are  adjusted  to  the  magnitude  of  the  amount  of  material  being  reacted  at  each  event 
When  less  than  1/F  mole  of  a  particular  species  is  left,  it  is  removed  from  the  RSV  and 
the  remaining  amount  of  moles  are  transferred  to  its  right  hand  side  neighbor  within  the 
RSV. 

3.3  Algorithm  and  software 

A  simplified  algorithm  of  the  numerical  code  developed  for  this  project  is  shown 
in  Figure  1.  The  complete  program  listing  with  variable  definitions  can  be  found  in 
Appendix  A.  Except  for  some  graphical  interface  functions  included  in  Macintosh 
specific  subroutines,  the  code  is  ANSI  FORTRAN77  standard.  The  implementation  of 
the  stochastic  logic  and  the  manipulation  of  reaction  vectors  are  carried  out  in  a  modular 
way  by  an  extensive  use  of  subroutines.  The  program  can  handle  a  reactor  space 
including  more  than  3000  different  species  (molecular  weight)  as  long  as  the  size  of  the 
RSV  remains  under  60  000  integer  elements.  These  limitations  result  from  memory 
resource  considerations  and  could  be  overcome  on  a  more  powerful  machine.  The  initial 
data  required  for  each  simulation  are  gathered  within  a  text  file  read  by  the  code  at  the 
beginning  of  the  execution.  The  template  for  this  data  file  is  included  in  Appendix  B.  The 
code  has  been  used  to  simulate  stochastic  processes  comprising  up  to  five  different 
events.  However,  more  complex  situations  could  also  be  studied  with  minor 
modifications.  The  results  of  the  simulations  are  saved  in  two  output  files:  poly.res  and 
polyds.res.  The  first  file  includes  the  time  evolution  of  Mn,  Mw  and  X  while  the  second 
file  reports  the  molecular  weight  distribution  as  a  function  of  the  conversion. 
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Assign  all  the 

Calculate  the 

code  control 
variables 

of  each  event 

* 

V 

Built  the  reaction 
vectors:  RSV(), 

lwvo.swvo 

..i - 

1 

Call 

Randomize() 

1 

/^np 


Input  initial 
data  on  the 
reactive  system 

A 

(  START  ) 


/ Loop  while  N. 
><  X<Xmax  > 


Print  RSV(), 
SWV(),LWV0 
to  file 


Select  the  reaction 
sites 


Perform  the  reaction) 
and  build  the  new 
molecule  vector 


Update  SWV() 
RSV()  and  LWV() 


■>:  n 


X,  iter,  N  species 

<--I  CaU  1 

Mw,Mn  to  file 

■  PLOT  MW  ■ 

and  screen^^— 

FIGURE  1  -  Simplified  Flowsheet  of  the  Computer  Code 


P154742.PDF  [Page:  29  of  78] 

UNCLASSIFIED 

17 

The  program  execution  begins  by  evaluating  an  initial  probability  factor  for  each  event 
before  entering  the  main  loop.  Within  this  loop,  successive  random  events  will  be  chosen 
and  performed  until  achievement  of  the  desired  level  of  conversion.  Whenever  the  amount 
left  of  a  particular  species  becomes  smaller  than  the  reaction  event  increment,  the 
subroutine  "REDUCE"  is  called  to  eliminate  this  species  from  each  polymerization 
vectors.  An  on-screen  progress  of  the  simulation  is  given  to  the  user  by  plots  of  Mn,  Mw 
and  MWD.  After  completion,  the  user  is  prompted  to  print  out  the  final  screen  layout.  An 
example  of  printed  screen  is  given  in  Appendix  C. 

One  of  the  options  given  to  the  user  by  the  code  control  variables  is  to  obtain  Mw  -  X 
and  Mn-  X  plots  for  the  considered  system  as  calculated  from  Macosko  formula  (Refs. 
20,21).  For  a  system  having  reactive  sites  of  the  same  reactivity,  these  curves  should 
coincide  with  the  same  plots  obtained  by  stochastic  simulations. 


4.0  SIMULATIONS 


4.1  Comparative  polyurethane  systems 

In  order  to  evaluate  the  efficiency  of  the  stochastic  algorithm  for  the  simulation  of 
polyurethane  polymerization,  three  different  systems  were  considered:  HTPB-TDI,  PPG- 
HDI  and  GAP-IPDI.  Their  chemical  structures  are  shown  in  Figure  2.  In  the  HTPB-TDI 
reaction,  the  two  isocyanate  groups  on  the  TDI  do  not  exhibit  the  same  reactivity.  Hence, 
three  stochastic  events  can  take  place:  1-  reaction  between  an  OH  group  and  a  NCO  para- 
substituted  (p-NCO);  2-  reaction  between  an  OH  group  and  a  NCO  ortho-substituted  (o- 
NCO)  or  3-  a  clock  event.  For  the  PPG-HDI  polyurethane,  NCO  at  both  ends  of  the 
hexane  chain  have  the  same  reactivity,  but  the  polypropylene  glycol  contains  primary  and 
secondary  OH.  Thus,  the  set  of  events  in  the  reactor  space  is  formed  by:  1-  reaction 
between  a  primary  OH  and  a  NCO  group;  2-  reaction  between  a  secondary  OH  and  a 
NCO  group;  or  3-  a  clock  event.  Finally,  the  GAP-IPDI  system  is  very  similar  to  the 
HTPB-TDI  reaction,  even  though  the  difference  in  reactivity  of  the  NCO  groups  is  much 
less  significant  than  with  TDI.  For  each  of  the  previous  polyurethane  systems,  the 
probability  of  events  was  estimated  from  kinetics  rate  constants  found  in  Table  in.  For 
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some  instances,  these  constants  were  not  available  for  the  conditions  of  the  simulations 
and  they  had  to  be  extrapolated  from  data  reported  in  the  references  of  Table  m. 


HO-(CH(CH  3)  -  CH2-  0)n-H 

PPG 


OCN-(CH2)-NCO 


HDI 


HO-(CH2  -CHCCH^MO-CHtCH^  )-Ciy-  OH 

GAP  +  1PDI  OCNH2C 


FIGURE  2  -  Initial  Components  for  each  Polyurethane  System 


TABLE  m 


Rate  Constants  used  for  Probability  Estimation 


reaction  1 

reaction  2 

reaction  type 

ki  (50  °C) 
(1/mol/min) 

reaction  type 

k2  (50  °C) 

(1  /mol/min) 

ref. 

HTPB-TDI 

p-NCO,  OH 

.142 

o-NCO,  OH 

.0035 

22 

PPG-HDI 

prim.  OH,NCO 

.0015 

sec.  OH,  NCO 

.0005 

23,24 

GAP-IPDI 

ring  NCO,  OH 

.0004 

meth.  NCO,  OH 

.0003 

23,24 
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4.2  Effect  of  polymerization  reaction  variables 

The  manipulated  variables  acting  on  the  polyurethane  formation  reaction  were,  for 
the  purpose  of  this  study,  the  temperature  and  the  NCO/OH  ratio.  It  was  expected  that  the 
temperature  would  affect  not  only  the  rate  of  reaction,  but  also  the  molecular  weight 
distributions,  for  systems  where  a  difference  existed  in  the  reactivity  of  reactive  groups 
for  the  polyol  or  the  polyfunctionnal  isocyanate.  On  the  other  hand,  the  stoichiometry  of 
the  reaction  determines  the  level  of  conversion  where  a  rapid  increase  of  the  molecular 
weight  is  observed.  For  each  of  the  systems  described  in  section  4.1,  a  simulation  was 
carried  out  for  three  levels  of  NCO/OH  ratio:  0.8, 1.0  and  1.2.  The  effect  of  temperature 
was  investigated  only  for  the  HTPB-TDI  reaction. 

The  results  of  the  simulations  are  represented  by  two  curves:  i)  a  plot  of  Mw/Mw° 
against  X  with  an  overlay  for  the  time  -  conversion  response,  and  ii)  the  final  molecular 
weight  distribution  at  conversion  level  of  80%.  These  curves  are  shown  in  Figs.  3  to  8 
for  each  of  the  three  polyurethane  systems. 


FIGURE  3  -  Effect  of  NCO/OH  ratio  on  HTPB-TDI  polymerization  at  50  °C 
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FIGURE  4 -MWD  for  HTPB-TDI  polymerization  (X=0.8)  at  50  °C,  NCO/OH=1.0 
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FIGURE  5  -  Effect  of  NCO/OH  ratio  on  PPG-HDI  polymerization  at  50  °C 
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FIGURE  6 -MWD  for  PPG-HDI  polymerization  (X=0.8)  at  50  °C,  NCO/OH=1.0 
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FIGURE  7  -  Effect  of  NCO/OH  ratio  on  GAP-IPDI  polymerization  at  50  °C 
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FIGURE  8  -  MWD  for  GAP-IPDI  polymerization  (X=0.8)  at  50  °C,  NCO/OH=1.0 


5.0  EXPERIMENTAL 

In  order  to  carry  out  the  stochastic  simulations  of  the  polyurethane  binder 
formation,  the  initial  molecular  weight  distribution  of  HTPB,  PPG  and  GAP  prepolymers 
had  to  be  determined  experimentally  by  gel  permeation  chromatography  (GPC)  based  on 
a  polystyrene  calibration.  The  polyols,  PPG  (Aldrich  Co.),  HTPB  R-45M  (Arco 
Chemical)  and  GAP  (3M)  were  used  as  received  from  the  supplier.  The  water  content  of 
these  materials  was  measured  by  Karl-Fisher  titration  and  it  was  found  to  be  below 
0.05%  for  each  polymer.  The  characterization  of  the  prepolymer  is  reported  in  Table  IV. 
For  the  PPG,  a  correlation  was  available  to  estimate  the  real  molecular  weight  from  the 
results  obtained  by  polystyrene  calibration  (Ref.  25). 
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TABLE  IV 

Characterization  of  Liquid  Prsp.Qlymgis 


Measured  by  GPC 

[polystyrene) 

Supplier  1 

Mn 

Mw 

polydispersity 

Mn 

Mw 

polydispersity 

HTPB 

6220 

12000 

1.93 

2600 

7000 

2.7 

PPG  1 

3100 

3200 

1.03 

3000 

3250 

1.1 

GAP 

2416 

3810 

1.57 

2000 

3200 

1.6 

1:  Corrected  by  correlation:  MppQ=5.731  Mp$0-714 


Moreover  a  limited  number  of  experimental  molecular  weight  determinations  were 
performed  on  HTPB-TDI  and  PPG-HDI  systems  to  assess  the  effectiveness  of  the 
computer  code  to  simulate  these  processes.  Since  both  systems  were  known  to  behave  as 
an  A2+  B2  reaction,  most  of  the  polymerization  products  were  soluble  in  an  organic 
solvent  at  high  conversion  levels.  The  isocyanates  were  obtained  from  Aldrich  Co.  The 
polymerization  was  carried  out  in  small  beakers  on  a  50-g  scale  within  a  temperature- 
controlled  oven.  During  the  reaction,  five  or  six  3-g  samples  were  removed  at  unequally- 
spaced  intervals  of  time  .  Each  sample  was  quenched  with  an  excess  of  diethylamine 
(DEA)  to  neutralize  the  remaining  isocyanate.  After  adding  the  DEA,  the  level  of 
conversion  in  the  sample  was  measured  by  13C-NMR  spectroscopy  (Ref.  22)  and  the 
MWD  was  estimated  by  GPC  in  4  successive  WATERS  Ultrastyragel®  columns  of 
grade  105  (stock  #  10574),  104  A  (stock  #  10573),  500  A  (stock  #  10571),  and  100  A 
(stock  #  10570). 


6.0  DISCUSSION 

At  this  point,  it  should  be  recalled  that  the  interest  of  performing  Monte-Carlo 
simulations  on  simple  polymerization  reactions  was  to  obtain  the  MWD  of  these  systems 
at  any  level  of  conversion  before  the  gel  point.  A  first  step  toward  the  validation  of  this 
work  regards  the  comparison  of  lower  moments  of  the  MWD  predicted  by  the  software 
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with  those  calculated  from  the  Macosko  and  Stockmayer  approach  (Refs.  18, 20  and  21). 
For  a  system  where  each  reaction  site  is  known  to  have  about  the  same  reactivity,  as  for 
GAP-IPDI,  it  was  found  that  stochastic  results  are  in  good  agreement  with  the  curves 
predicted  by  the  more  conventional  approach. 

It  should  also  be  noted  that  unequal  reactivity  will  only  affect  Mw  since  Mn 
depends  only  on  the  number  of  molecules  reacted,  not  their  length.  An  example  of  this 
effect  is  shown  in  Figure  9  for  stochiometric  HTPB-TDI  polymerization  at  30  °C.  In  this 
case,  the  NCO  site  in  para  position  is  expected  to  be  50  times  more  reactive  than  the  ortho 
substituted  NCO  (Ref.  22).  Thus,  there  is  a  chain  elongation  effect  taking  place  that 
delays  the  weight  average  molecular  growth.  In  accordance  with  this  hypothesis,  the 
results  of  the  simulation  are  perfectly  coherent 


FIGURE  9  -  Effect  of  asymmetrical  reactivity  on  HTPB-TDI  polymerization 
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FIGURE  10  -  Comparison  between  measured  and  calculated  degree  of  polymerization  for 

HTPB-TDI  reaction  at  NCO/OH  =1.0 


FIGURE  1 1  -  Comparison  between  measured  and  calculated  degree  of  polymerization  for 

PPG-HDI  reaction  at  NCO/OH=1.2 
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As  discussed  in  section  5,  experiments  have  been  carried  out  in  the  laboratory  to 
compare  the  simulation  results  with  real  polymerization  products.  This  comparison  was 
based  on  Mn  and  M\y  for  HTPB-TDI  and  PPG-HDI  systems.  Figures  10  and  1 1  show 
that  there  exists  a  fairly  good  agreement  for  both  systems  between  calculated  and 
measured  molecular  weights  except  for  high  levels  of  conversion.  The  differences 
observed  may  be  associated  with  the  fact  that  the  measurements  were  obtained  from  a 
polystyrene-calibrated  GPC,  since  the  factor  between  actual  and  measured  molecular 
weight  may  change  for  higher  molecular  weight.  This  is  also  suggested  by  the  fact  that 
this  difference  is  more  important  for  the  PPG-HDI  curve,  where  the  non-linear  correlation 
available  was  developped  for  low-molecular-weight  polymers. 


FIGURE  12  -  Comparison  between  measured  and  calculated  MWD  for  PPG-HDI 

system  at  NCO/OH=0.8 
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Based  on  the  fact  that  the  calculated  Mn  and  Mw  appeared  to  satisfy  experimental 
evidences,  it  is  justified  to  assume  that  the  whole  MWD  is  reasonably  well  simulated  by 
the  software.  However,  one  of  the  problems  with  the  manipulation  of  the  calculated 
MWD  is  to  reconstruct  a  continuous  distribution  from  a  series  of  discrete  mass  fraction 
and  individual  molecular  weights.  This  can  be  accomplished  by  usual  regularization  code 
with  equation  [12]  as  a  constraint  Even  without  an  extensive  mathematical  treatment  a 
simple  comparison  between  calculated  and  measured  MWD,  as  shown  in  Figure  12, 
reveals  an  acceptable  similitude. 

Another  important  issue  to  ensure  a  meaningful  simulation  was  the  frequency 
factor  determining  the  amount  of  moles  reacting  during  a  single  event.  Based  on  an  initial 
quantity  of  100  moles  of  material,  it  was  found  that  F  should  at  least  be  greater  than  100. 
This  implies  that  10000  reaction  events  were  needed  to  acheive  complete  conversion. 
Because  a  greater  F  value  translates  into  more  species  within  the  RSV,  the  calculated 
MWD  is  therefore  closer  to  the  real  physical  system  where  an  extremely  large  number  of 
different  polymer  chains  are  found.  If  only  the  Mn  and  Mw  values  are  of  some  interest,  a 
smaller  F  will  still  give  results  expected  from  Stockmayer's  equation  with  a  smaller 
computational  effort.  A  typical  run  time  for  a  simulation  was  about  30  minutes  on  a 
68040  based  personal  computer. 


7.0  CONCLUSIONS 

Stochastic  modeling  has  been  used  to  simulate  the  polymerization  of  a  diol  with  a 
diisocyanate  for  polyurethane  systems  similar  to  those  found  in  solid  rockets  motors. 
The  computer  code  developed  for  this  purpose  has  proved  to  be  reliable  for  the  prediction 
of  the  MWD  as  a  function  of  the  level  of  conversion  (advance  of  the  reaction).  This 
information  about  the  molecular  structure  of  a  reacting  system  can  be  used  to  estimate  the 
rheological  behavior  of  these  materials.  This  would  allow  a  more  detailed  optimization  of 
propellant  formulations  as  well  as  a  better  design  capability  for  their  related  industrial 
processing. 

However,  to  achieve  truly  useful  results,  the  same  kind  of  algorithm  would  have 
to  be  adapted  for  systems  with  significant  cross  linking,  that  is,  systems  with  overall 
equivalent  functionality  greater  than  2  for  the  polyfunctional  isocyanate  and  polyol.  This 
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would  require  minor  modifications  in  the  software  for  the  encoding  mode  of  information 
about  the  molecular  structure  of  species  in  the  reaction  space. 

In  summary,  it  has  been  demonstrated  that  Monte-Carlo  simulations  offer  a 
powerful,  yet  simple,  mean  of  investigating  the  development  of  advanced  polyurethane- 
based  energetic  materials  with  a  minimum  of  experimental  information.  Accordingly, 
there  is  a  great  interest  in  pursuing  the  development  of  a  numerical  code  that  will  allow  the 
simulation  of  more  complex  systems,  more  representative  of  in-service  propellant 
formulations.  The  work  described  in  this  memorandum  will  be  followed  by  the 
experimental  rheology  characterization  of  solid  rockets  propellants  and  the  modeling  of 
the  behavior  of  these  systems  will  be  accomplished  using  the  results  presented  here. 
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ptri  •  FrontRindomO 

CALL  ECRAK2,2*.5«.229.*59.''iPtr2-,rind0l"J5 

CALL  POLICEC21,*,*) 

froieWtVENEICNT 

CALL  (toveToCVol2(5),Vol2C2»» 

CALL  Oro.Te«tCfro»eS,Vat2C»).VAL2aS)) 
f  rcxeS-  -  ' 

CALL  *bv«ToCVol2CS5,Vo\2C3S)) 

CAL L  Dr o*Text(f roteS , VoI2C»),VALZC25» 
froteS-'X*<  -  ' 

CALL  MoveTo(ValZ(S),VoT.2(5R)) 

CALL  Oro«T*xtCf roteS, Vol2(e).VAL2C2S)> 
fraseS-'X  •  * 

CALL  *>veToCVal2CS),Vol2C65)) 

CALL  OraoTextCf roteS ,Vol2C») .VAL2C2®)) 
froteS-’naifcre  total  d  etwees:' 

CALL  H3veTeCVol2C5),V°T2C»e)) 

CALL  Ora«TextCf roseS, VoT2C»),VAL2C2S5) 
froteWtoille  du  vecteur  RSV:' 

CALL  ttoveTo(Vol2(S)  ,Val2C95)) 

CALL  DroeTextCfroteS ,Vol2C») ,VAL2(2S)> 

CALL  WTTTIXECptrl) 
tenpss-'eee 

CALL  SHOWTTH£(t*«iptS,ptrl) 

CALL  ECRANC2, 2M,  1,479, 319,  ’  *  ,ptrZ , m ndo«Z) 

CALL  :NITPLOTCptr2,l,2Se,319,A79,0.,l.,l»ee.,5.1EA>l, 

♦  SX2,$Y2,L2,H2) 
froteS-*Mn  or  to  vt  X* 

CALL  ttovtToCVoT2(12«),VAL2Cie)> 

CALL  POlICE(21,12,1) 

CALL  ForeCoLorCVoU<4»9)5 

CA.L  Dro«Text(frotfS,Val2C*),VAL2C2)) 

CALL  ForeColorCVal*C33)) 

CALL  Dro*Text(froseS ,Vol2CZ),VAL2C3)) 

CALL  ForeColor<Val*C341)) 

CALL  Dro*Text(frateS ,Vol2(5),VAL2(4)) 

CA.L  ForeCclorCVol<C33)) 

CALL  DroeTextCf roseS, Vo12C3),VAL2C5)) 

CALL  ECRANC2, 2S#,32», *79,639,  "  ,ptr3,m'idw3) 

CALL  DiITPL0TCptr3,32e,2Se,639,«79,ie.,3.ES,e.,ie..«,S)C3.Sv3, 

*  13, H3) 

CALL  POLICE C21, 12,1) 
froieS-"t*<  Otitributson" 

CALL  kt>vtToCVol2C12e),VAL2Cie» 

CA^L  Dro*Text(f roseS, Val2(0) ,VAl2(16)) 

CALL  Select*; ndow(VAL*CPtrC)) 

CALL  SttPortCV*L*Cptre» 

CA.L  ImtCurtorO 

KRITEC6,*)  'APPUYEZ  SUP  RETOUR  POUR  LIRE  UN  FICHIER  DE  DOSSEES" 
PA,SE 
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C  Main  program  starts  here 


*f>0LY-5«. 

TIME-4. 

X-«. 

C  . . . . •• . •••••••••••••••#•••••••••••< 

c  LECTURE  DCS  DONNEES  A  PARTIR  D’lW  FICHIER  TEXTE 


OPENClf.FILE-** ,STATUS-*0U>’ .FORM-’ UNFORMATTED’) 

READClf,*)  tit\«S 

READU®  *)  F 

READ<1®,*)  XMAX 

READC1®,*)  NCOOH 

READU®,*)  NISO 

READU®,*)  MOH 

READU®,*)  PROP 

READU®.  *>  DT 

READU®  ,*)  « 

READU®,*)  K2 
READ(1®,*)  K3 
READU®,*)  K4 
READU®,*)  DWIIME 
READU®,  *>  MAC 
READU®, *)  NT 
RUDU®,*)  MONO 
READU®,  *)  IBODY 
DO  1-1 ,  NT 

RE ADC 14,*)  POPULC2,I),SWVCI) 

END  DO 
CLOSE (14) 

C  . . . . . . . . . 

C  OVERTURE  D’UN  FICHIER  DE  RtSULTATS 

C  . . . . . . . . . 

UNIT-11 

OPE  NC 11 ,  Fit  F- '  poly .  res  *,  STATUS- •  OLD  *,  FORM- ’  FORMATE  D ' ) 
RRITECUSIT,*)  title! 

WRITE (UNIT, •)  -DISTRIBUTION  INITIALE' 

WRITECUhIT,*)  ’  t  «  ,Qt*' 

DO  1-1 , NT 

WRITE (UNIT ,1®)  I,POPULC2,I),SWV(I) 

END  DO 

OPE  N(  12 ,  F IL  E- ’  pol  yds .  res ST  ATUS- ’  OLD  *,  FORM- ’  FORMATE  D  ‘ ) 
WRITEC12 ,*)  title! 

CALL  MDLECULEOiT  ,MDNO,IBOOY  ,MR,MN,SN,  .FALSE  .  ,6) 

CAU  PLOTDS(NTf SX3,SY3,L3,H3,ptr3) 


C  . . . 

c  FORMATION  DES  VECTEURS 

c  . . . . . 

C  fV*CTIONAL  croups  table 

c 

C  t  FUNCTIONAL  CROUP  PRIM.  0« 

C  1  FUNCTIONS  GROUP  SEC.  OH 

C  2  FUNCTIONAL  GROUP  NCO 

C  3  FUNCTIONAL  GROJP  URETHANE 

C  4  MONOMER 

C  5  ISOCYANATE  BODY 

C  6  FUNCTIONAL  GROUP  NCO' 

C  . . . . 

C  F OWING  the  REACTOR  SPACE  VECTOR 

C  POLYMER  CHAINS  FIRST 

miso-apoly*  ncooh 
fl-0 

00  I-l.NT 
Fl-FL+I 

IFCNOH.EQ.l)  THEN 
ASVCFL>1 
ELSE 

rsvcfl>-« 

END  IF 
FL-FL-1 

R5V(F L)-4*POPULC2 . I)/MONO 

Fl-Fl^l 

RSV(FL>1 

P0PULC2 ,  D-WULC2 ,1)+34 
END  DO 


C 


ISOCYANATECS)  FOLLOWS) 
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Fl-FL+l 

RSV(FL>2 

FL-FL+1 

RSV(FL>-5 

FL-FL41 

RSVCFD-2 

IF  (NIS0.EQ.2)  THEN 
FL-FL+1 
RSVCFl>2 
FL-FL+l 
RSVCFL>5 
FL-Fl+l 
RSVCFL>€ 

EM)  IF 

C  FORMATION  OF  THE  SPECIES  WEIGHT  VECTOR 

C  SUR  UNE  BASE  DC  MPOIY  DE  POLYMER  0>0  HI  SO  MOLES 

C  D’YSOCTANATE  C*CO) 

C  . . ••••-••••••••••••*••• . . 

c  HYDROXYL  CHAINS  FIRST 
DO  1-1,  NT 

SWVG  )-SWVCI) /SN»»*>OL  Y 

END  DO 

C  ISOCYANATE  (S)  FOLIOW(S) 


NT-NT+1 

SWV(NT).*ISO*PROP 
POPUt  (2 ,  NT>IB00Y+*4 
IF  (NISO.EQ.2)  THEN 
•fT-MT+l 

SrvCKT>MISO*a .  -PROP) 

POPULC2,WT>IBOOY^84 

EM)  IF 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


FORMATION  OF  THE 
UNA  TAKE 
HNK(l)  -  Ct-O 

LlMC(2>  (1-4) 
lINAC3)  - 

LIKA(4>  (2-5) 
LlNK(S)-  (5-6) 
lINK(6)-  C3-4) 
LlNK(7).  C3-5) 


LINK  WEIGHT  VECTOR 

PRIM.  OH-MONOMER 
SEC.  Oh- MONOMER 
MONOMER -MONCVER 
NCO- ISOCYANATE  BODY 
NCO' -ISOCYANATE  BODY 
URETHANE -MONOMER 
URETHANE -ISOCYANATE  BODY 


DO  1-1,7 
LINK(X)-*. 

EM)  DO 

SPECIE-1 

LAST-RSV(l) 

DO  1-2, Fl 
NCM-RSVCI) 

I F ((LAST . EQ . f . AND . NOW . EQ . 4) . OR . ( LAST . EQ . 4 . AND . NOW .  E  Q .  6)  )  THEN 
LINKCD-LINKCD-fSWYCSPECIE) 

EM)  IF 

IFCCLAST.EQ.B.AND.NCm.GT.6).0R.CLAST.GT.6.AND.NCm.EQ.e))  THEN 
LlHKCl>LlNKCl>SlfYCSPE  CIE> 


IF((LAST.EQ.l.AND.NOW.EQ.4).OR.CLAST.EQ.4.AND.NOW.EQ.l))  THEN 
LIXK(2)-LINKC2>S*V(SP£CI£) 

EM)  IF 

IF(CUST.EQ.l.AND.N0W.GT.6).0R.(LAST.CT.6.AND.N0*.EQ.l))  THEN 
LINKC2>LINKC2>SWVCSPECIE) 

EM)  IF 

IF  (N0W.GT.6)  LIMCC3>aMKC3>SWVCSPECIE)»CFLOATCNOW)/4.-10 
IF  CO-AST. EQ. 2  .AND.  NOW.EQ.S)  .OR.  (LAST.EQ.5  .AND.  N0W.EQ.2)) 
-»  LINNC4)-LINKC4)-»SWVCSPECIE) 

IF  ((lAST. EQ. 6  .AM).  NOW.EQ.S)  .OR.  (LAST.EQ.5  .AND.  N0W.EQ.6)) 
4  LIKKCS)-LlNKC5)4SFrYCSPECIE) 

IF  (CLAST. EQ. 3  .AND.  MOW. IQ. 4)  .OR.  (LAST.EQ.4  .AND.  N0W.EQ.3)) 
4  LINKC6)-LINKC6)4SWVCSP£CIE) 

IF  (CLAST. EQ. 3  .AND.  NOW , GT .6)  .OR.  (LAST>€  .AND.  NOW.EQ.3)) 

4  LINKC6)-LINKC6)4SWVCSPECIE) 

IF  ((LAST. EQ. 3  .AM).  MOW.EQ.5)  .OR.  (LAST.EQ.S  .AND.  NOW.EQ.3)) 
4  llNKC7)-lINK(7)4S»fY(SPECIE) 

if  ((last. eq. 4  .am>.  now.eq.d  .or.  clast.eq.s  .and.  now.eq.2) 

♦  .OR.  CLAST.EQ.S  .ANO.  M0W.EQ.6))  SPECIE-SPCCIE41 
IF  ((LAST. GT. 5)  .AND.  (NOW.EQ.l))  SPECIE-SPECIE41 
LAST-MOW 
EM)  X 


WRITE (UNIT,*)  ‘REACTOR  SPACE  VECTOR:" 
WRITECUNIT,11)  (RSVCD.I-l.FL) 
WR1TECUNIT,*)  ‘SPECIE  WEIGHT  VECTOR:" 
WRITE  (UNIT,  12)  (SWCD.I-l.NT) 
WRITECUNIT,*)  *LI>*  WEIGHT  VECTOR" 
WRITE (UNIT, 13)  CLlNK(I), 1-1,7) 
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IF  Q04AX  GT.NC00HM96.)  WAX^COOH*100. 

itenwx-INTCXMWF) 

CALL  MDLECULECWT.M^.IMOY.fcW.W.SN,  .FALSE  .  .NISO) 

IF  (MAC)  THEN 

CAL  L  MACOSKOC I»OOY  ,*« ,  NCOOh  ,  L2 ,  H2 ,  SX2 .  SYZ , pt r2 ,  NT ,  NI SO .  XMAX) 
E®  IF 

CALL  MDLECULECNT.FCW.IEOOY.Ffc.MN.SN,  .TRUE.  ,NIS0) 

OtDXCl)-INTC.l*LZ-X*SX2>l 

OLDXCZ>-OU)xa) 

OU>Y  ( 1)-H2- INT( .  1*H2^«*SY2) -1 
OLDY(2>H2-IN7(  .1*H24***SY2)-1 

CALL  PL0T>«CX>^>SX2,SY2.0LDXa),OLDYa),l2,H2,ptr2IA«91Z) 

CAL  L  PLOT**  (X  ,  SX2 ,  SY2 ,  OLOXC2  ) ,  Ou D Y  C2 ) ,  L2 ,  H2 ,  pt r2 , 341 , 2  ) 

WITECUNIV)  •increment  -  *,1/F,’  Holes’ 

WRITE(UNIT,*)  ’NT  FUG  ITER  X  Hn  M*  time  (mi n)* 
WRITEC6,*)  'increment  -  \1/F,'  moles’ 

CALL  RANOOMI2EO 

CALL  POP(FL,WT,MONO,IROOV) 

1-1 

DO  9HILECI.LE.NT) 

IF(5WY(I) .LT.l/F)  THEN 

CALL  REDUCE (FL, NT ,1 .SITE, LINK, 1) 

ELSE 
1-1*1 
END  IF 
END  DO 

PClOCK-l./DT 

iter«l 

DO  WHILE (i ter. LE.itermax) 

C  SELECT  AN  EVENT 
SITE -9 

DO  WHILE  (SITE . EQ  6) 

P1-K1*L1NX(2)*l:nXC4)»F/1S«  . 

F2-*2»LINX(2)»LINXC5)-F/I5e . 

P3-X3*LlNXa)»LINKC4)*F/lSe . 

P4»X4,l1NK(1)*lINK(5)*F/150 . 

PTDT-P14P24P34P-UPCL0CK 

P5-PUOCK/PTOT 

P4_P4^TOT 

P3-F3/PTOT 

P2-P2/PTOT 

Pl^l/PTOT 

CHOIX^NDC) 

IF  CCH0IX.GT.CP1*P24P^P4))  THEN 
TIME -TINE -DT 

ELSE  IF(CH0IX.CT.P14.P2-P3)  THEN 
SITE-6 

o^e 

ELSE  IF CCHOIX , GT . Pl-PZ)  THEN 
SITE-2 

Qf-,  g 

ELSE  IF(CHOIX.G’.Pl)  THEN 

site-g 

OH-1 

ELSE 

SITE-2 

OH-1 
END  IF 

C  IF  C(LINX(4) . LE . 1/F) . AND . (SITE . E Q . 2))  SITE-0 

C  IF  ((LlNKCS).tE.l/F).AND.(SITE.EQ.6>)  SITE-0 

END  DO 

C  LOCATION  OF  THE  REACTION  SITES 

OECK-6 

DO  WHILE  (CHECK. EQ. 9) 

SPCO**« 

chpnt-lx>*ci*oh)*Rndo 

SPCNCO-0 

IF  (SITE .EQ.  2)  THEN 
NCOPNT-LINK(4)*RND() 

ELSE 

NCOPNT-LI*§C(5)*RND() 

END  IF 
SOMOH-0. 

somnco-9. 

C  RECHERCHE  DE  L  ESPECE  OH 

IF  CRSV(l).EQ.OH)  SOMOH-SOMOi*SW^(l) 

IF  (SOMOH.GE.OHPNT)  THEN 
flT>-l 
SPCOfvl 

DEIUT-POPUL(l.l) 

HPACE--1 
END  IF 

IF  (RSV(POPULCl.l))  EQ.OH  .AND.  SPCOH.EQ.6) 
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♦  SOMOH-SOMOK*S*VC1) 

IF  (SOCH.CE.OHP^)  THEN 

ftn-POPULCl»l) 

SPCOH-1 
DEBUT-1 
*>ACE-1 
END  IF 
1-2 

DO  WHILE  ((SPCOH.EQ.t).ANO.I.lE.NT) 

IF  (RSV(POPULO,I-1)+1).EQ.OH) 

*  SOOH-SOOH^SWVCI) 

IF  CCSOM0H-C«>NT)/F.CT.-l.E-3)  then 
f\n-POPUl(l,I-l)+l 
SPCCH^I 

DEIUT-P0FULC1,SPC0H) 

»ACE— 1 
END  IF 

IF  CRSVCP0PUIC1,D).EQ.0H)  THEN 
SCMOH-SOMOW+SWCI) 

IF  C(SOMOH-OHPNT)/F.GT.-l.E-3)  THEN 
fan-POPULCl.I) 

SPCOi-I 

DEBUT -POPULCl ,SPC0H-1)^1 
HPACE-1 
END  IF 
END  IF 
1-1*1 
END  DO 

C  RECHERCHE  DE  L  ESPECE  NCO 

IF  CRSVCD.EQ.SITE)  SOMNCO-SOMNCOfSWVCD 
IF  CSOfcttCO.GE.NCOPNT)  THEN 
LDEBUT-1 
SPCNCO-1 
L FIN-POP UL Cl,  D 
PACE-1 
END  IF 

if  ccrsvcpopjici,D).eq.site).and.cspcnco.eq.0))  THEN 
SONCO-  SO**CO*SNV(l) 

END  IF 

IF  CSOMNCO.GE.NCOPNT)  THEN 

ldebut-popjlci.D 

SPCNCO-1 
LFIN-1 
PACE— 1 
END  IF 

1-2 

X  WHILE  CCSPCNCO.EQ.0)  -AND.  I.LE.NT) 
if  crsvcpopulci,i-D*D.eq.s:te) 

.  S0MNC0-SO4NCOSWVCI) 

IF  CCSOMNCO-NCOPNT)/F.GT.-l.E-3)  THEN 
LDEIUT-POPULCl,I*l)-*-l 
spcncg-i 

LF IN-POP UL ci, SPCNCO) 

PACE-1 
end  if 

IF  CRSV(POPULa,I)).EQ.SITE)  THEN 
SO**  CO-SOMNC  0+  SWVC I) 

IF  CCSCKNC0*NC0PNT)/F.GT.-l.£-3)  then 
LD£IUT4>0PUlCl,I) 

SPCNCO-I 

L  F IN-P0PULC1 , SPCNCO- 1 )+l 
PACE— 1 
END  IF 
END  IF 
1-1*1 
END  X 
CHECK-1 

IF  CSPCOH.IE .•  .OR.  SPCNCO. LE.0)  THEN 
CHECK-« 

END  IF 

IFCCSPCOH.  EQ .SPCNCO) . AND .CSW(SPCOH) . LT . 2/F))  THEN 
CHECK-0 
END  IF 

END  X 

C  ICI  ON  FORME  LA  NOUVELLE  MOLECULE  DE  SA  PARTIE  Oh 
COOT-1 

X  I-OEBUT,fin-HPACE ,HPACE 
N£*VCCXNT)-RSVCI) 

COUNT-COUNT^l 

END  DO 

SOH-COUNT 

NE*TVCC0UN7>3 

C  PARTIE  NCO 


CXNT-COUNT^l 
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00  I-LDEBUT+PACE ,  IFIN,  PACE 
NEWV(COUNT)-RSV(I) 

COUNT  -COUNT  ♦  1 
END  DO 

COUNT  ^OONT-l 

C  PERFORMING  THE  REACTION  AND  ADJUSTING  VECTORS 

5INCSPC0H)-SWVC$PC0H)-1/F 
IF  CSWCSPC0H).LT.1/F)  THEN 

IF  CSPCOH.LT. SPCNCO)  5PCNC(VSPCNC0-1 
CALL  REDUCECFLjNT fSPCOH, SITE , LINK, OH) 

END  IF 

SfW(SPCNC0)-S1WCSPCNC0)-l/F 
IF  (SWVCSPCNC0).LT.1/F)  THEN 

CALL  REDUCECFL, NT, SPCNCO, SITE, LINK, OH) 

END  IF 

HNKCl4OK>-LINKClf0H>-l/F 

LINK(6)-LINK(6>1/F 
LINK(7)-LINKC7>1/F 
IF  (SITE.EQ.2)  THEN 

LINK(4>LINKC4)-1/F 

ELSE 

LINKC5)-LINKC5)-1/F 
END  IF 


C  CA.CUL  DU  POIDS  MOLECULAIRE  OE  LA  NOUVELLE  ESPECE 

poiDs-e 

DO  1-1, COUNT 

IF  (NEW(I).EQ.4  .OR.  NEWVCI) . GT .6)Th£K 
POIDS-MONO-NENVCD/A+POIDS 
END  IF 

Ic  CNENV(I).EQ.S)  POIDS-IBOOY-POIDS 
IF  CNEWV'CD.EQ.OH)  P0IDS-17+P0IDS 
IF  CNE*N(I).EQ.SITE)  POIDS-42-POIDS 
IF  CNEWCI).£Q.3)  P0IDS-59+P0IDS 
END  DO 

C  ON  CHERCHE  L*EXISTENCE  DC  CE  POIDS  MOLECULAIRE  DANS  POPJL 

INSERT-© 

JUMLAU-© 

1-1 

DO  IMILE  CI.LT.NT) 

IF  (CP0IDS.EQ.POPUlC2,I)).AND.CJlWEAU.EQ.e))  JUMEAU-'I 

IF  CP0:DS.LT.P0PULC2,I))  then 
INSERT-I 
I-AT 
END  IF 
I-I-rl 
END  DO 

IF  CINSERT.EQ.©)  INSERT-NT+I-NISO 
CHECK-© 


c  s:  GROOVE 

IF  CJUMEAU.NE.6)  THEN 
DO  JJ-JUMEA^. INSERT 

IFCCCCNEIlVCl) .  EQ.RSVCPOP  JLC1 ,  JJ-1>1)) .  AND .  CNEtfV(COJNT) .  EQ. 

*  RSVCPOPoLCl.JJJJJJ.OR.CCNE^CU.EQ  RSVCPOPJLCl,JJ))).AND.(NEifVCCO 

•  UNT) . EQ.RSV(P0PULC1, JJ-lJ+l)})) .AND.CCnECK.EQ.©))  ThEN 

CHECK-1 

SWVCJJ)-SWV(JJ>1/F 
END  IF 
Ett  DO 
END  IF 

C  NON  TROUVE 

IF  (CHECK. £Q.«)  THEN 

DO  I-FL,P0PU-C1,INSERT-1>1,-1 

ASVa+COUNT>RSVa) 

END  DO 

DO  1-1, COUNT ,1 

RSVCP0PJLC1, INSERT- 1)+I)-NEWVCI) 

END  DO 

00  IWfT.INSERT.-l 
5WV(I*1>-SKV(I) 

POPULC2,I-l)-POPULC2,I) 

POPULCl  ,1+D-POPULCl  ,I)+COUNT 
END  DO 

.  SNV(INSERT)»1/F 

P0PJL(2 ,INSERT)-P0IDS 
POP  JLC1 ,  INSERT)-P0PUL  Cl  .INSERT  -1>C0UN7 
NT-NT-t-1 
FL-FL+COUNT 
END  IF 

XCNKC6)/W>OLY/2 . 

CALL  MDlECULECNT.MDNO.IBOOV.MR.MN.SN,  .TRUE. , NISO) 

C  »r;t*CUhlT,»)  SN.Se^SC.'NCOCH-FLQATCiteO/F.iter 
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C  FIN  t)E  I'EVENEKENT 

CALL  S«tPortCVAL4(ptrl)) 

CALI  EFFACE (7®, 13, 3®, It) 
froseS-ISTRGSCiter) 

CALL  MDveToCValZC7e),Val2C2e)) 

CAL L  OronTextCf rose! ,Vfl\2C®) , VALZC1S)) 
CALL  EFFACE (105, 73, 30, 10) 

frose$-lSTRGSCNT) 

CALL  *tov«ToCVol2Cie5),Vol2CM)) 

CAL L  Dro*T«xt Cf rose* , Vol 2C0) , VALZ(15)) 
CALL  EFFACE  (100,8*.  »,1®> 

froseS-lSTUGJCFL) 

CALL  MoveToCVol2C10e), Vol 2(95)) 

CALL  DrowTextCf roseS , VolZ(0) , VALZ(IS)) 


IF  (H0DCX,8.82).LT.l/F/iee)  THEN 
CALL  EFFACEC35.2M3.1®) 
froseS-RSTRGS(MK) 

CA.L  HDveToCVol2C35),Val2C35)> 

CALL  Dr«mText(f roseS, Val2(0),VAL2(10)) 

CALL  EFFACEC35.43.43.ie) 

frt>ScS-RSTRG$0*O 

CALL  *>veToCVol2(35),VclZC5e)) 

CALL  DrcmTextCfraseS ,Val2C0) , VAL2C10)) 

CALL  EFFACEC35.58.43, 10) 
froseS-RSTRGS(X) 

CALL  *>veTo(Val2(35),Val2C6S)) 

CA.L  Dro»T«xtCfrosel .Vol2C0) .VAL2C5)) 

CA^  t  PlOTV*(X  ,MA ,  5X2 ,  SY2 , 0lDX(1  ) ,  0LD?a  ) ,  L2 , H2 .  pt rZ  ,409 , 2  ) 
CA_  L  PL0T>*f(X  ,Mh , SX2 , 5Y2 , OIDXC2 ) . OwDY (2  ) ,  L2 ,  h2 ,  ot  r2 , 34 1 , 2) 
NRITE(UNIT,14)  NT.FL.iter.X.MN.MA.TIME 
hejre-lNTCAlNTCTI(«  /60 . )) 
froseS-ISTRGS(hejre) 

\f (neure.lt .10C)  f roseS- *@‘//fraseS 
\f (he jre.lt .10)  froseS-'® ’//f roseS 
ternpsS-froseS(l:3)//* :  * 
hfrjre- JMOOCI KT(A1  NT  CTIME ) ) , 60) 
froseS-lSTRGSCheure) 
ifCheure.lt. 30)  fraseS-'0’//f roseS 
tempsS-teffpsSCl :  4)//f  roseS(l  :2) 
heure-INTCCTIMc  -  AINT(nME))*6e . ) 
froseS-ISTRGS(heure) 

*f(heure . It . 10)  f roseS-’® *//f rose! 
tempsS-teepsSCl  :6)//' :  V/fraseS 
CALL  SrtOTTlMECtefnpsS.ptrl) 

EM)  IF 

IF  CMOD(X,0.1).IT.  l/F/100)  THEN 

CA.L  PtOTDSCNT,SX3,SY3,L3,H3,ptr3) 

IF(IMPRIME)  THEN 

C  NR ITEC6.*)  “ESPECE  0H-,CNEWYCI),I-1.S0H-1),1 

C  NRITE(6,«)  -ESPECE  NCO"  .SHE  ,(N£WV(I),3-S0h>»1,C0lNt) 

C  NRJTE(6 ,*)  -KEN  MJ.ECULE\(NENV(I),I-1  .COUNT) 

C  WRITEC6.15)  1,P0PULC2,1)*SWVC1),CRSVCI),I-1,P0PJLC1.D) 

C  J-2 

C  DO  NHILECJ.LE.NT) 

C  NRITEC6.15)  JjPOPUlC2,J),SIWCJ),CRSVCI).I*POPJLC1,J'1>1.POP 

C  *  ULCl.J)) 

C  J-J-l 

C  END  DO 

NRITEC12,*)  "NT-  -.NT,"  FL-  ",FL,"  X-  *,X,"  ITER-  ",\ter 
DO  J-l.NT 

NRITEC12.16)  J,POPULC2,J),S«VCJ) 

END  DO 
END  IF 
EM)  IF 
tter-iter+l 
END  DO 

CALL  PLOTDSCNT,SX3,SY3,L3,H3 ,ptr3) 

WRITE(6,*)  ‘TERMlNE ‘ 

NRITECUNIT,*)  ’DISTRIBUTION  FINALE* 

NRITECUNIT ,15)  1 . POP JL(2 , 1) , SNVC1), CRSV(I) , 1-1 . P0PUL(1 ,1)) 

J-2 

DO  NHILE(J.LE.NT) 

WRITEOAIT ,  15)  J  .POPUL (2 , J) , SNV(J) , CRSV(I)  ,I-P0PUL(1 , J-1>1 , 
*•  POPUL Cl. J)) 

J-J+l 
EM)  DO 

WRITECUNIT,*)  "LINK  NEIGHT  VECTOR" 

WRITS (UNIT, 13)  CLINK(I), 1-1,7) 

CuOSECH) 

CL0SEC12) 

WRITE (fe,‘)  ' VOULEZ-VOUS  INFRIKR  ?  (OUI-l.NQN-e)  ’ 

READC5.-)  I**5 
IFCIMP.EQ.l)  THEN 

CALL  GPOuwpCVol4Cptrl),Vol4Cptr2),Vol4Cptr3),0) 

END  IF 
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PAUSE 

1»  FORMAT  ('  \7X,  IS,',  MS,',  \  F6.3) 

11  FORMAT  (’  •,7*CK,1X)> 

■2  FORMAT  (•  MKF7.J,  ZX)) 

FORMAT  f  \1KF7.2,  ZX)) 

4  FORMAT  C  \I4,  V.IS.V.IS,  VF4.3,  V.F6.»,  V,F6.6,\\F6.1) 

15  FORMAT  c  \ZX,IS,4X,I6,4X,  F6. 3, 24(14, IX)) 

16  FORMAT  C  ' ,ZX,li,*X,l6,*X,  F6.3) 

EM> 


C 


*•*  FIN  DU  PWXRAA44E  PRINCIPAL 
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C  S0U5-PR0«AmE  MOOIFIAMT  LA  POLICE  B'IMPRESSION 

. . . . 

SUBROUTINE  P0LICECN1,N2,N3) 
integer**  N1,K2,N3 
CALL  TextFontCVol2CNl)) 

CALL  T«xtS\*«CVol2(NZ)) 

CALL  T«xtFaceCVd2(N3)) 

E® 


C**** . 

C  5tr255  function 

C  M OKt  0  Posed  string  from  o  FORTRAN  string 

CHARACTER*25€  FUNCTION  Str2$5(string) 
CHARACTER*^)  string 

Str255  -  CHARCtENCTRIK*tring)))//string 

RETURN 

E® 
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C  . . ........ 

c  SOUS  -PROGRAM  CREEANT  UNE  FENFTRE 

SUBROUTINE  ECRAN(ren, top, left, bottom, right .title, ptr, window) 

RECORD  /Rect/  wrect  !  rectangle  defining  window 

RECORD  /»indo«*ecord/  window  !  storage  for  window 

INTEGER*-*  top, left, bottom, right, ptr, nm 

CHARACTER***  title  !  title  string  for  windows 

CHARACT£R*256  5tr2S5  !  fortron  to  pascal  string  function 

EXTERNAL  Str255  !  defined  ot  end  of  this  file 

•rect .top  -  top  1  where  to  put  second  window 

wrect. left  -  left 

wrect .bottom  •  bottom 

wrect. right  •  right 

ptr  •  NewRindo»»Cwindow,wrect,Str255Ctitle),vallC-true.), 

♦  vol2Crw),voUC0)lvollC.folse.),vol4Ce)) 

CALL  Select«indowCVAL4Cptr)) 

CALL  SetPortCVAl4CptO) 

RETURN 

£® 


3/23/95  16:34  AM  Hardisk: FORTRAN  STUFF :polymer.f .§  Page  1 


c  SOUS -PROGRAMME  CONVERTISSANT  UN  ENTIER  EN  CHAINE  DE  CARACTERES 

c  ••••• . . . 

FUNCTION  ISTRGSCW>«RE) 

INTEGER**  M0M6RE ,L 

CHARACTER*15  ISTRGS 

CALL  N^ToStringCVALACWOMBRE), ISTRGS) 

If  (NCNBRE .NE .6)  THEN 

l-intcalogikabscfloatcno»«re))))+2 

ELSE 

L-2 

E»IF 

IF  CNOMBRE.LT. 6)  THEN 
L-Ul 
EMDIF 

ISTRGS-ISTRG$a:L) 

RETURN 

EW 

C  SOUS-PROGRA*€  CONVERTISSANT  UN  REEL  EN  CHAINE  DE  CARACTERES 

C  ••••••»■ . 

FUNCTION  RSTRGS(NOMBRE) 

CHARACTER*3€  RSTRGS, RESTS 
INTEGER**  ENTIER, L 
REAL**  NOMBRE 
ENTIER-INTCAINTCNCWBRE)) 

CALL  »Aer.ToStringCVAL*CENTIER),RSTRGS) 

IF  CNCMBRE.NE.6.)  THEN 

L-INTCAL0G1«CABS(N0MBRE))>2 

ELSE 

1-2 

E®IF 

RSTRGS-TRIMCADJUSTL(RSTRGS)) 

ENn  E  R-ABSCINTCCNOMB  RE  -  AINTCMOWRE))  *1666)) 

CAlI  hLm»TeStrir>gCVAL*CEKTlER)  .RESTS) 

DO  I-INT(AL0G16CFL0AT(ENTI£R)>1),2 
RESTS-*6“//R£STS 
E®  DO 

RSTRGS-RSTRGSC2 :  L)//CK1RC*6)//RESTSC2  :■ *) 

RETURN 
E  » 
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C  CE  SOUS- PROGRAM  «NERE  LE  FORMAT  D’UN  GRAPH 
C 

C  *  -  UNDO*  NlACER 
r  XI-  UPPER  left  corner  pixel  x 
Yl-  UPPER  left  comer  pixel  v 
v  X2-  BOTTOM  ri girt  corner  pixel  x 
C  Y2-  IOTTOM  right  corner  pixel  y 
C  XHIh,  XMAX,  YKIN,  YHAX  -  scole  definition 
C  GRID  -  t:  No  grid,  GRID-1,  flrid  lines 

c  . •••••••••••••• . •••••••••••••••••••••••••••••••••••••••• 

SUBROUTINE  INnPLOTCptr.Xl.Yl.XZ.YZ.XMIN.JOUX.YMIN.YMAX.G.SX.SY, 

*  larfleuTjhouteur) 

INTEGERS  ptr,Xl,Yl,X2,Y2,G(lorgeur,Houteur,ox,o> ,XOlD,YOlD 
REALM  XXIX , XMAX , YMIN , YMAX , PROPX , PROPY , SX , SY 
CHARACTER*3®  NS 

CHARACTER*2S6  ISTRGS , RSTRGS  !  integer/REAL  to  string!  function 
CHARACTERISE  Str255  !  fortran  to  pascol  string  function 

EXTERNAL  Str2S5, ISTRGS, RSTRGS  !  defined  ot  end  of  this  file 
RECORD/Rect/r 

CALI  SetPortCVAUCptr)) 

lergeur-X2-Xl 

houteur-Y2-Yl 

PROPX-.* 

PROPY-.* 

SX-PROPX  •  lor  geur /CXKAX  -  XWN) 

SY-P ROPY  •  hout  eur/CYMAX  -  YMIN ) 
ox-HlNTCCl>PR0PX)/2*lorgeur) 
oy-*JINTCCl-PR0PY)/2»houteur) 

XOLD-ox 

YOlD-oy 

r.top^jy 
r. left-ox 

r  .  bottom-oy+NINT  CPROPY*ho  uteur) 
r .  right-ox-NINTCPROPX*lorfleur) 

CALL  PenSizeCVAL2C2),VAL2C2)) 

CALL  FraneRectCO 

CALL  PenSireCVAL2(l),VALZCl)) 

CALL  POLlCEC21,9,e) 

CALL  NaveTo(Vol2(ox-25),Vol2Coy+2)) 

NS- ISTRGS  Cl  NTCYMAX)) 

CX.L  OropTextCNS, Vol2CB) > VAL2CINT(ALOG10CTMAX))+1)) 

CALL  MoveToCVol2Cox-5>,Vol2Coy*NINT(PROFY*nojtejr>10)) 

NS-ISTRGSCINTCXMIN)) 

CA.L  Oro*Text(NS,Vol2(e),VAL2CINTCAL(X;i0CXV:N))Tl)) 

CALL  McveToO'elZ(ox-25),Vol2(oy-*-NINTCPROPY»nQjtcjr>2)) 
ns-istrgscintcyvin)) 

CALL  0ro*Text(NSlVol2Ce),VAL2(lNTCAL0GieCYMIN))45)) 

CALL  ktoveToCVoi2(ox»IWT(PROPX»lorflejr>?0), 

-  VaUCoy4NlffT(PROPY*nouteur)4l0)) 

NS-ISTRGSCINT(XHAX)) 

CAlL  Dro*.Text(NS,Vol2Ce),VAL2CIKTCALOG10CXVtAX)).l)) 

IFCGEQ  1)  THEN 
1-2 

DO  fP-ILECI .  LE  .8) 

CA.L  MoveToCVol2(ox),Vol2(py4NINT(PROPY*houteur/20“I))) 

CA.L  LineCVol2ClMCPROPX*\argeur)-l),Vol2C0)) 

CALL  MoveToCVol2Cox4NINTCPROPX*lorgeur/10-:)),Val2(oy)) 

CALL  L\neCVol2C0),Vol2CINT(PROPY-nouteur)-l)) 

1-142 
END  DO 
E *j  IF 

CALL  MoveToCVol2Cr. left), VolZCr. bottom)) 

RETURN 

END 
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C  CE  SOUS -PROGRAAWE  CALCUL  1*00,  **CN)  ET  U  QT£  TOTALS  (MOLES) 
C 

C  INPUT:  WT,»CNO,I*OOY  OUTPUT:  »*,**, SN 


SUBROUTINE  ICLECULECWT.MONO.IBOOY.MK^SN, CHOIX, NISO) 
logical  CHOIX 

vnteger*2  RSV(W©ee),NEWVC2tee0) 
irtteger*4  POPUtCZ.ieee) 
integer**  NT,*ONO,IIOOY,FIN,NISO 
reol*4  SirVC3We),*r,HN,SNfP0I0S 
reol*l  StM.SWC 
CtMON  POPUL.RSV.NEPTV.SWV 
SW4-B. 

SN-4. 

StiC  ‘ft. 

IF  (CHOIX)  THEN 

F IN-NT 

ELSE 

FIM-NT-NISO 
END  IF 
DO  2-1, FIN 

POIDS-*lOAT(POPUl(2 , I)) 

S>*uPOIDS*SttV(I>St* 

SN-S*fVCI)4SN 

S*C-P0IDS*P0IDS*S*VCI>SW2 

END  DO 

**-*EA^CSWC/SNM) 

**-REAc(S*<)/SN 

RETURN 

END 
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C  CE  SOUS 'PROGRAMME  CALCUL  WE  HiTRICE  COMPRENANT  LE  POIDS 
C  MOLECJLAIRE  D€  CHACUNE  DES  ESPECES  ET  LEUR  POSITION  TERMlNA. 
C  DANS  LE  REACTOR  SPACE  VECTOR 

C  . . . 


SUBROUTINE  POP(FUG ,NT,MONO,IBOOY) 
integer*2  RSV(6B00e)  ,NEWV(2W00) 
trrteger*4  P0PUL(2 ,3*») .SPECIE ,m, FLAG 
integers  NT.kCNO.IBODY.lAST.NO* 
real *4  SWV(SWe) 

COAON  POPUL.RSV.NERV.SWV 
LAST-RSV(l) 

m-* 


SPECIE-4 
DO  1-1, FLAG 
NOB-RSV(I) 

IF  CM0N.EQ.4  .OR.  N0N.GT.6)  fc*-MN*MONO*NO*/4 
IF  CN0R.EQ.5)  **«4**IBOOY 
IF  (MOB.EQ.l.OR.  NOR  .EQ.B)  MB-174** 

IF  ((NOP  EQ.2)  OR .(NO*. EQ.6))  HK-42+m 
IF  (N0R.EQ.3)  **-59*>* 

IF  ((LAST. £0.4  .AND.  NOW.LE.1).OR.(LAST.EQ.5  .AND.  N0*.EQ.2).0R. 
(LAST.GT.6  .AND.  NOft.LE ,1)  .OR.  (LAST .EQ.5  .AND.  NW.EQ.6))  THEN 
SPECIE-SPECIE*! 

P0PUL(1,SPECIE)-I 
P0PUK2,  SPECIE)-** 

END  IF 
LAST-NO* 

END  DO 
KT-SPECIE 


RETURN 

END 
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C  C£  SOUS -PROGRAMME  ELDCNE  LIS  ESPECES  DONT  LA  QT|  PRESENTE 
C  EST  INFtRIEURE  A  l*  1/F 


SUiROUTINE  REDUCE (F LAG , NT , I , SITE ,  LINK ,OH) 
integer-2  R5VC«e«e)  ,N£WC2ie«) 

integer-4  POPUIC2 ,3000) , NT , FLAG , SHIFT , START , I , FOH , FHCO , SIT E , OH 
rm\* 4  ST/C 3*00), AMOUNT, LINK(Z) 

CO440N  POPUl,RSV,NEWV,SWV 
FOH-e 
FNCO-O 

AMOUNT-SRVCI) 

SHin-POPULCi.i) 

START-1 

IF  CX.GT.l)  THEN 

SHI  FT-POPUL  Cl ,  I) -WUl  Cl  ,1  - 1 ) 

STARTO*OPULC1,I‘1>1 
END  IF 

IF  CRSVCSTART).LI.l)  FOH-1 
IF  CRSVCPOPUlCl,I)).LE.l)  FOH-FOH^l 
IF  CRSVCSTART).EQ.SHE)  FNCO-1 
IF  CRSVCP0PULC1. 1)5. EQ. SITE)  FNCO-FNCO+1 
DO  J-I.NT-l 
SWV(J>SWCJ-1) 

POPJL  Cl , )  >POPUL  Cl ,  J^l) -SHI  FT 
P0PULC2 ,1>P0PUIC2,  J+l> 

END  DO 

DO  J-ST ART, FUG-SHIFT 
RSVCJ)-RSVCJ*SHIFT) 

END  DO 

FUG-FLAG -SHI  FT 
NT-WT-1 

IF  CRSVCSTART).LE.l)  foh-foh-i 
IF  CRSVCPOPUlU,D).LE.l)  FOH-FOH-1 

IF  crsvcstartj.eq.she)  fnco-fnco-i 
IF  CRSVCP0PULC1,I)).EQ.SITE)  fnco-fnco-i 
SWVa>SNVCI>-  AMOUNT 

lin>;ci-*Oh>linkci*oh)-CR£AlCFOh))#amdjnt 
IF  CSITE.EQ.Z)  THEN 

LINKC4)-LINKC4)-CREALCFNC0))*AM0UN't 

ELSE 

LINKCS)-LINKC5)-CREALCFNC0))*AMOUKT 
END  IF 

RETURN 

END 
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c  . . ***** 

C  C£  SOUS -PROGRAMME  TRACE  LA  DISTRIBUTION  DE  POIDS  MOL  ECUL AIRE 

C  . . . .......MM . I . . 


SUBROUTINE  PLOTDSCNT.SX.SY.L.H.ptr) 
integer *2  RSVC60ee0),NEirV(20ee0) 
integer**  POPULC2.3®®®) ,ox,oy,IX,IY 
integer**  NT,L,H,ptr 

reel**  SINC3®«)  .SX.SY.LEHAX.closseCSB)  ,cot 

RECORD /Rect/r 

COWON  POPUl,RSV,«KV,SWV 

o*-INTC.l*l>l 

oy-INTC.l*H>l 

r.top-oy+2 

r.left-ox+5 

r.bott»INTC.9*H)-Z 

r.right-INTC.9*L)-5 

CALL  SetPortCVAL*Cptr)) 

CALL  ErcseftectCO 
cot-* .*77/5® . 

DO  J-1,5® 
cla*seC3>-®. 

END  DO 
DO  1*1, MT 

J-INTCCALOC10CFLOATCPOPUI(2,I))>1.>/Cot) 

IFCCJ .GT .0). AND. O .LE  .56))  closseCJ>closseCJ>S#VCI) 
E*>  DO 
IEHAX-®. 

00  1-1,5® 

IF  (classeCD-GT.lEMAX)  LEHAX-closseCI) 

END  DO 
DO  1-1,5® 

lX-©x+INTCCFLOATCI)*COt)*.S*FL0AT(iy*.*77) 

CALL  MoveToCVoI2CIX),Val2(IN7C.9*H)-l)) 
IY-INTCclosseCI)/LEMAX*9.*SY) 

IF  CIT.U-CINT (.**«)))  THEN 

CALL  L\neTo(VoI2CIX) ,Val2CINT( .9*H)-I v-1)) 
END  IF 

E  DO 
RETURN 
END 


c  . . ••'**•*•* . •*•••• . ••••••' 

C£  SOUS -PROGRAMME  TRACE  LE  POIDS  MOtECULAIRt 


SUBROUTINE  Pt0TV*f(X,MK,SX,SY,0LDX,0LDY,L ,H,Dtr,cojlejr,grCS) 
integer**  ox,oy,lX,IY 

integer**  L ,H,ptr ,OLDX,OLDY,couleur,gros 

r#ol*4  X,HI,SX,SY 

Cx-lNT( .l*l)+l 

Oy-INTC . 1*H)+1 

CALL  Set  Port  CVAL*Cptr)) 

lX^x-INT(X*SX) 

IY-*.-o>-INTO*rSY) 

IF  (iY.GE.oy)  THEN 

CALL  ForeColorCValACcouleur)) 

CALL  PenS;xeCVAw2(gros),VAL2Cgros)) 

CALL  MoveToCVoTZCOLDX) ,Vol2(OLDv)) 

CALL  LintToCVoI2CIX),Vol2(IY)) 

CALL  ForeColorCVflL*(3®)) 

E*>  IF 

OlDX-IX 

OLDY-IY 

CALL  PenSizeCVAL2Cl),VAl2(l)) 

RETURN 

END 
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C  S0US-PR06RAWE  SERVANT  A  ALIM  ENTER  LE  GENERATED  OE 
C  NOMBAE  AlEATOIRE  AVEC  LA  DATE  ET  L'HEURE  DU  SYSTEME 

C  . . . 

SUBROUTINE  RAJOQMIZEQ 
INTI  GERM  GETAS 
INLINE  CGETAS-z^eBd') 

IntegerM  sec 

RECORD  /QmckdrowGlobals/  QOGlobols 
POINTER  (p.QDGlobols,  QOGlobols) 

P-QOGlobols  -  L0NGCGETA5Q)  -  SIZEOFC/QinckdrawGlobols/) 
CALL  GctDoteT\meC*ec) 

QOGlobols .rondSeed  -  sec 

RETURN 

END 

C  SOUS -PROGRAM  FONCTION  RETOURNANT  UN  NQM8RE  ALEATOIRE 

FACTION  UNDO 
REALM  AND 

R«)-C32767 .  -FLOAT  CRondomO))/65534 . 

RETURN 

EH) 


. . . . . . ... 

C  SOL'S -PROGRAM  EFFACANT  UNL  CHAMP  DANS  UNE  FENETRE 

...... ........ . ****** . . . . . 

SUBROUTINE  EFFACECX.Y.L.H) 

integerM  X,Y,L,H 

CALL  ForeColorCVol4(3e)> 

CALL  PenS\zeCVAL2(6),VAL2C9)) 

CALL  »toveToCVol2CX)fVolZCY)) 

CALL  L\neTo(Vol2CX.L) ,ValZCY)) 

CALL  ForeColorCVo\4(33)) 

RETURN 

END 
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C' 

c 

C' 


SOUS -PROGRAMME  ESTIHANT  Nn  et  M*  selon  Mocosko . 

SUBROUTINE  MACDSkOCIBOOY,*Ba,NCOOH,1.21H2,SXZ,SY2,ptr2,NT,NIS01E) 

integers  Rsvc»«ee),NEirYC2BBee) 

integer**  P0PUlC2,i*®e). NT, WSO, BORNE 

integer**  IBOOY , 01DX( 2) , 0LDY£2) , 12 , H2 ,pt r2 

r«ol*4  IB.  ,*B»I ,  NCOOH ,  X ,  SX2 ,  SY2 ,  FE ,  GE ,  BODY ,  R 

reel**  MAFAF,»efBr,SAF,SBF1iiia,**>,mop,i*p,MAF2AF,M8F28F 

real**  4feW,  CM ,  SKV(3*i»; ,  E 
COBCN  POPUl,RSV,NENV,S«V 

MAFAF-B. 

MAFZAF-B. 

SAF-B. 

DO  I-1.NT-NIS0 

HAFAF  -MA'AF  .POPUU2  ,I)*SWVCI) 

MAF2AF  J«F2AF^>0PUL  (2 ,  l)*POF>ULC2 ,  I>*S*fva) 

SAF-SAF+SEFYCI) 

END  DO 
WFBF-B. 

WFZIf-B. 

SBF-B. 

DO  I-NT-NTSO*l,NT 

WFIF-MBFgt  J>OPUlC2,I)*SWVCI) 
!«F2tFJ«t2BF.P0F’UlC2.I)*POPUlC2,I)*SNVa) 

SBF.SBF.SrYCI) 

END  DO 
FE-2. 


GE-2. 

■o-MAFAF/SAF/FE 

mb-»«FBF/SBF/GE 

map  JAAF2AF/SAF/FE 

■PP-NBF2BF/S6F/GE 

0*A-FE*«s 

»«-GE*irt 

BOCY-FlOATCIBOOY) 

R-1/NCOOh 

t*^(HAFAFJ«FBF)/CSAF«SBF) 

MB-(Cl/NCOC")***II*>*iI-BODY*BODY5/Ck«:/NCOOr.»BOOY) 

0LEXC1>TNTC.1*L2>1 

OLDXC2>-OLDXCD 

OLDY(l)-H2- INTC .  1  *H2»***SY2)  - 1 

OLDY  C23-H2  -  INTC .  1  *H2»*«,  *5  Y2  )  - 1 

BORNE-INTCE) 

DO  I»l,BORNE,l 
X-FlOATCIVIBB. 

(*,-OULFAF*l«FBFVCSAF«SBF-FE*X*SAF) 

CA'.L  PlOT)*»CX,F»<,SX2,SY2,O.DXC2),OtDYC2)1L21H2,ptr2,3B,l) 
m.fR*X*map*X*mbc>/CR*X*«w»X*'”P) 

HU»l.X*X*R*CX*CFE-10*9«,0«*X*R*CGE-l-)*W(A!A-2*GMA*GNB)/ 
CX*R*mo-X*mbVa.  -X*X*R*CFE-1  ,)*(GE  -1  )) 

CA„L  Pt0Tk*lCX,mi,SX2,SY2,0lDXCl),0:D>Cl),F2,H2,ptr2,33,l) 

END  DO 


RE'uRN 

END 


P154742.PDF  [Page:  63  of  78] 


UNCLASSIFIED 


3/23/95  10:4«  AH  Hardisk : FORTRAN  STUFF: polymer. f .§ 


C  SOUS -PROGRAM  INITIAWT  L'HORLOGE  VIRTUEL 

c . . . . . 

SUBROUTINE  INlTTIHE(ptr) 
integer**  ptr 
RECORD/Rect/r 
r. top-15® 
r. bottom-186 
r. left-15 
r .right-115 

CALL  SctPort(VAL4Cptr)) 

CALL  ForeColorCVal*C33)) 

CALL  PenSxreCVAL2C2),VAL2C2» 

CALL  FrameftectCO 
CALI  ForeColorCVoHC33)) 
r .top-133 
r  .bottom-156 
r. left-15 
r.  right -115 
CALL  FraneRectCO 
r. top-153 
r  .botto"^177 
r. left-18 
r  .nght-112 

CALL  PeaS\2eCVAL2Cl),VAL2Cl)) 

CALL  FrameRectCO 

CALL  MoveToCVol2Cl8),Vol2Cl46)) 

CALL  POl2CEC21>12,6) 

CALL  OrcmTextC* VIRTUAL  CLOCK*  tVol2(0),VAL2C13)) 
CALL  P0LICEC21,9,6) 

RETURN 

END 


c— . . . . . 

C  THIS  SUBROUTINE  DISPlAVS  ™e  VIRTUAL  TIME  INTO  NlNDO*  Ctrl 

. . . . . 

SUBROUTINE  SHOBTIMECTEHPSS.ptrl) 
integer**  ptrl 

chorocttr*9  TEMPSS 
RECORD/Reet/r 
r .top-15* 
r. bottom-176 
r. left-19 
r.  ngbt-111 

CALL  SetPortCVAL*Cptrl)) 

CALL  EroseRectCr) 

CA.L  H3veTo(Vol2(26)IVol2(172)> 

CA^L  POLlCEC21,18 ,1) 

CALL  ForeColorCVoi*C265)) 

CALL  DrcmTe*taEvrSS,Vol2C6),VAL2C95) 

CALL  POlICEC21,9,«> 

RETURN 

Etf> 
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C  This  Subroutine  prints  (litnop)  winbons  1  ,2  ond  3 

SLBROUTIRE  GPDurc>CgrafPortl  ,grofPort2  .grafPortS ,  bitShope) 
IMPLICIT  WOKE 
INTEGER*!  errCode 

IKTEGER*<  grafPortl,  bitShope,  bMopPtr,  rectPtr 
IKTEGER#4  grofPort2 , grofPort3 
errCode  -  • 

CALL  PrOrvrOpenO 
errCode  -  PrErrorO 
IF  (errCode  .RE.  »)  GOTO  188 
CALL  PrCtlCcll(VAL2CiPrDtvCtl>, 

A  VAL4ClPrReset),  VAL4(»),  VAL4(8)  5 

errCode  -  PrErrorO 
IF  CerrCode  -RE.  8)  GOTO  188 
bttopPtr  -  LOCCgrofPortl)  r  2  !  bitMop  offset 

rectPtr  -  LOC (grofPortl)  +  16  !  portRect  offset 

CALL  PrCtlCoUCVAL2CtPr»ltsCtI), 
i  VAL*(bMopPtr),  VAL4(rtctPtr),  VAL4(bitShape)  ) 

bMopPtr  -  L0CCgrofPort2)  ♦  2  !  bitMop  offset 

rectPtr  .  LOC(grofPort2)  *  16  !  portRect  offset 

CALL  PrCtKoUCVAL2CiPrBitsCtl), 

A  VAL4(bMopPtr),  VAL4(rectPtr),  VAL4(bitShope)  ) 

bMopPtr  -  L0CCgrofPort3)  ■»  2  !  bitMop  offset 

rectPtr  -  L0CCgrofPort3)  «  16  !  portRect  offset 

CALL  PrCtlCall(VAL2(lPrBitsCtl), 
i  VAL4(b*topPtr),  VAL4(rectPtr)t  VAL4(bitShape)  ) 

errCode  -  PrErrorO 
IF  (errCode  .RE.  8)  GOTO  180 
CALL  PrCtUoUCVAL2(\Pr0evCtO, 

A  VAt4(LPrPogeEnd),  VAL4(8),  VAL4(e)  5 

errCooe  »  PrErrorQ 
IF  (errCooe  .RE.  8)  GOTO  188 
CA.L  PrCtlCoUCVAL2(iPrDev<tl>, 

A  VAL*aPrDocClose),  VAL4C6),  VAL*(8)  ) 

errCode  -  PrErrorO 
IF  (errCooe  .RE.  8)  GOTO  188 
CA-L  PrDrvrCloseO 
errCode  -  PrErrorO 
188  contirue 
RETURR 
END 
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Template  for  input  data  file 
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Titre  de  la  simulation 

P  (tel  que  masse  1  4v4nement=l/F) 

conversion  totaled 

ratio  NCOOH 

nombre  d' isocyanate 

nombre  de  site  OH 

prop,  de  1 ' isocyanate  1 

incr4ment  en  temps  (min) 

R4activit4  du  site  NCO  #l-OH  sec. 
R4activit4  du  site  NCO  #2-0H  sec. 
R4activit4  du  site  NCO  #l-OH  prim. 
R4activit4  du  site  NCO  #2-OH  prim. 
Reduction  du  RSV() 

Impression  i  1‘dcran 
Prediction  Macosko 
nombre  de  classe 
MW  du  monomfere 
Mw  de  1' isocyanate 
GPC:  Mw,  qt4  en  mole 


Chalne  de  caractdres*72 
r4el 
r4eL 
r4el 
entier 
entier 
r4el 
r4el 
10. 

1. 

10. 

1. 

.TRUE.  or  .FALSE. 

■TRUE,  or  .FALSE. 

•TRUE,  or  .FALSE. 

entier 

entier 

entier 

entier,  r4el 


I 
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Example  of  Screen  Output 
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£YENEMENT  *  18000 
m=  14501.414 
Ivtst  =  7724.297 
Xr  0.880 

nombna  total  despeces:  1297 
taille  du  vecteur  RSY :  256 1 9 


VIRTUAL  CLOCK 

1147:56:00 


Mw  Distribution 
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